DEPARTMENT OF OCEAN ENGINEERING 

MASSACHUSETTS INSTITUTE OF TECHNOLOGY 
CAMBRIDGE, MASSACHUSETTS 02139 



ROBUST RATE CONTROL SYSTEM DESIGNS 
FOR A SUBMERSIBLE 

By 

LAWRENCE JOHN DREHER 



O.E. 

S.M. (ME) 



COURSE 1 3 A , 2 
MAY 1984 





ROBUST RATE CONTROL SYSTEM DESIGNS FOR A SUBMERSIBLE 



by 



Lawrence John Dreher 



B.S., Rensselaer Polytechnic Institute 
(1977) 



SUBMITTED TO THE DEPARTMENTS OF OCEAN ENGINEERING 
AND MECHANICAL ENGINEERING IN PARTIAL FULFILLMENT 
OF THE REQUIREMENTS FOR THE DEGREE OF 
OCEAN ENGINEER AND MASTER OF SCIENCE IN 



MECHANICAL ENGINEERING 



at the 



MASSACHUSETTS INSTITUTE OF TECHNOLOGY 



May 1984 



© Lawrence » ler , 1984 




TV s 

£ ")3 








ROBUST RATE CONTROL SYSTEM DESIGNS 
FOR A SUBMERSIBLE 

by 

Lawrence J. Dreher 

Submitted to the Departments of Ocean Engineering and Mechanical 
Engineering in the Partial Fulfillment of the Requirements for the De- 
gree of Ocean Engineer and Master of Science in Mechanical Engineering. 



ABSTRACT 

Two robust rate control system designs are carried out for a sub- 
mersible (modeled by the NSRDC 2510 equations) in a turn using the Linear 
Quadratic Gaussian approach with L oop T ransfer Re covery. Separate com- 
mand channels allow the submersible to maneuver independently in the 
horizontal and vertical planes; the vehicle operator controls the bearing 
rate and depth rate through a joystick-like device. The fin conf iguration 
is the conventional cruciform stern without differential control. 

The first compensator design, called the r - 6 controller in this 
thesis, directly controls two vehicle state variables; pitch (0) and the 
vehicle angular velocity, r, about the z-axis. The other system, called 
the ip - z controller, controls yaw (or heading) rate (ip) and depth rate 
(z) directly. However, this design relies on linearized equations of yaw 
and depth rate to be employed by the compensator for state reconstruction. 

A tool for Kalman Filter loop shaping is developed in which state 
variables are scaled to provide good loop shapes and then recovered to 
get controllers that are robust. Both controllers are compared on the 
basis of performance in a nonlinear simulation. A robustness comparison 
is also conducted. 

Based on limited simulation data, this thesis concludes that the 
^ - z controller provides better control of depth rate than the r - 0 
controller. Bearing rate performance is essentially equal in both de- 
signs. However, the ip - z controller appears less robust in certain 
frequency ranges. 
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CHAPTER 1 



INTRODUCTION 



1 . 1 Background 

The submarine operator’s objective is to maneuver his vehicle 
freely in the ocean environment with respect to depth, heading, and 
speed. In some situations, such as transits and partrols, the operator's 
task is simplified to maintaining the submarine's depth, heading, and 
speed essentially constant. Such maneuvering requirements are easily 
accomplished by the crewmembers of the submersible . On the other hand, 
docking, turning, and performing evasive maneuvers to avoid ' navigation 
hazards or enemy torpedos are instances where the submarine's depth, 
heading, and speed may have to be varied suddenly. Obviously these 
maneuvers must be carried out safely and effectively. ‘However, there are 
many environmental and design obstacles that inhibit the operator from 
fully controlling this process. These impediments include a cruciform 
stern without differential control; cross-coupling of the planes of motion 
(due primarily to the large sail area) , variable hydrodynamic forces 
(vorticity, turbulence) , and a limited ability to "see" the environment 
(no windows) . 

Largely through experience and familiarity with the vehicle 
dynamic response, the conning officer learns that certain rudder angle/pitch 
angle orders will result in a certain heading rate and depth rate for a 
given speed. Rules of thumb have been developed for ordering the appro- 
priate amount of pitch angle to accomplish the required depth change and 
for an appropriate amount of stern plane angle to offset the vehicle's 



tendency to dive in a turn. It is particularly difficult for the conning 
officer to simultaneously command the turn rate and pitch angle due to 
the interaction of complex dynamics. The operator's control task be- 
comes increasingly difficult as the severity of the maneuver increases. 
One would therefore seek to eliminate much of the guesswork associated 
with performing these maneuvers by providing more substantial means of 
controlling the transient dynamics of the submarine. 

It is the intent of this thesis to demonstrate a procedure for a 
truly multivariable control system design that would provide the operator 
with a means of achieving desired submarine motion in depth and heading 
through rate control. The multiple input - multiple output (MIMO) de- 
sign methodologies based on the MIMO-LQG formulation of Stein and Doyle 
[1] will form the theoretical basis of this thesis. Specifically the 
novel MIMO-LQG/LTR approach proposed by Athans [2] for model based com- 
pensators will be followed in this controller design. Importantly, 
this methodology permits the design of dynamic compensators without the 
need for full state feedback. This will allow the development of a 
"robust" controller that is tolerant of modeling errors, nonlinearities, 
and noise. 

1 . 2 Outline of Thesis 

Chapter 2 first introduces submarine dynamics and deve lopes the 
linear models that will be used for two different control system designs 
at 30 knots. The last section of the chapter will provide the desired 
controller design specifications to be adhered to in the remainder of 
this design. 

In Chapter 3 the theory that provides the basis for the LQG-LTR 
process is introduced. 

Chapter 4 is the design procedure chapter for the two controller 
designs. A method to produce consistent Kalman Filter loop shapes is 
developed. 



12 



In Chapter 5 the obtained designs are evaluated through computer 



simulation and comparisons are made. 

Chapter 6 contains the summary and suggestions for further study. 
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CHAPTER 2 



THE NONLINEAR SUBMARINE MODEL 



2 . 1 Introduction to Submersible Modeling 



Submersible dynamics of motion and attitude can be described in a 



variety of coordinate systems. For the purpose of developing a model of 
submersible motion, the equations of motion are expressed in the body 
fixed axis because hydrodynamic forces and inertia are most readily com- 
puted in a ship reference frame. From the aspect of ship guidance and 
control, on the other hand, it may be desirable to describe motion, such 
as vehicle course and depth, in terms of an earth reference frame. 



tion of ocean vehicle motion. These represent equalities of Newtonian 
force and moment expressions, on the left hand side, and the so called 
dynamical response terms on the right hand side: 



General dynamical equations have been developed for the descrip- 




f (dynamical response terms) 



( 2 - 1 ) 



The general form of the force expression (Newton) is: 



— (momentum) 
at 
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where 



U 



G 



iu + jv + kw 



and the moment expression is 




— (angular momentum) 
at 



d 

dt 




G 



where 




(The definitions of u, v, w, p, q, and r are found in Table 2-1). 

The subscript G indicates the origin to coincide with a body-fixed 
coordinate system located at the vehicle center of gravity. However, the 
reference point is seldom taken there since the center of gravity moves 
with shifting weights inside the vehicle. Instead, the coordinate sys- 
tem is generally taken about the submersible center of buoyancy, and so 
provides a useful location for hydrodynamic estimates. This point is a 
function of vehicle geometry and is therefore fixed. The left hand side 
equations become somewhat more involved due to this coordinate system 
transformation and will not be discussed further. Details of these 
terms and simpl if cations commonly used for ocean vehicle dynamics can 
be found in reference [3] . 

The dynamical response terms of the right hand side of Eq. (2-1) 
are presented in Appendix A. These terms express the external forces 
and moments exerted on the vehicle by hydrodynamic, control surface, 
propulsion, and other effects. 

The force and moment equalities of Eq. (2-1) describe the six pos- 
sible degrees of freedom of the submersible. The three forces are in 
the axial, lateral, and normal directions, which give rise to motions 
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Table 2-1. Definitions of submarine states and controls. 



Submarine States 



State 


Definition 




Units 


u = 


-f 

Component of U in direction of 


x axis 


(ft/s) 


v = 


— b 

Component of U in direction of 


y axis 


(ft/s) 


w = 


->■ 

Component of U in direction of 


z axis 


(ft/s) 


P = 


Component of ft about x axis 




(rad/s) 


q = 


Component of ft about y axis 




(rad/s) 


r = 


Component of ft about z axis 




(rad/s) 


<f> = 


Angular rotation (roll) about the 
x axis 


(radians) 


e = 


Angular rotation (pitch) about 
y axis 


the 


(radians) 



Submarine Controls 

6r: deflection of rudder 

6s: deflection of stern plane 

6b: deflection of sail planes 

RPS : shaft revolutions per second 
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of surge, sway, and heave respectively. The three moment equations pro- 
duce moments and motion of roll, pitch, and yaw. Figure 2-1 shows the 
positive directions of forces, moments, motions, and control surface 
deflections. 

2.2 The Nonlinear Computer Model 

The nonlinear model used for this design was derived from the 
original NSRDC 2510 document "Standard Equations of Motion for Submarine 
Simulation" [4] . These equations have since been improved to include 
crossflow drag and vortex contributions. This model is installed on the 
computer system at the Charles Stark Draper Laboratories. The nonlinear 
model used in this thesis consists of 8 differential equations to describe 
the submarine dynamics. The six equations derived from the force and 
moment equalities account for the states u, v, w, p, q, and r. The de- 
pendence of these states on hydrostatic restoring forces about the pitch 
and roll axis and their kinematic relations result in two additional 
states, $ and 9. The definitions of these states are listed in Table 2-1. 

At this point, it will suffice to state that ship motion in open 
water is not sensitive to heading angle; on the other hand, the dynamics 
of a submersible are affected by changes in buoyant forces brought about 
by moderate depth excursions. These forces were not accounted for in 
the submarine model and, hence, depth is not included as a state in the 
system equations. By choice, propulsion plant dynamics were excluded 
in this model; the resulting propeller dynamics in an actual submarine 
vary, depending on the operating procedure of engineering plant personnel 
and would be difficult to model. Moreover, the model in its current ver- 
sion does not include either the actuator dynamics, or the actual angle 
rate limits of the control surfaces. Due to the importance of actuator 
dynamics on the control of actual submersibles , they will be considered 
as part of the augmented plant dynamics for this design in Section 4.1. 
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M,q 




STARBOARD VIEW 



\ 





Figure 2.1. 



Conventions showing positive directions of deflections, forces, 
velocities, and moments. 
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The resulting eight nonlinear equations as used in the computer 
simulation with the definitions of the hydrodynamic coefficients are given 
in Appendix A. These coefficients are approximated by direct measurement 
on a full scale or model submarine, or analytically. The accuracy of 
the dynamic response of the model is governed to a large extent by the 
accuracy of these hydrodynamic coefficients. Additionally, some forces 
have yet to be modeled mathematically, and other effects such as vortex 
shedding and separation effects are not possible to include in a linear 
model and are most likely the weakest point in the design model. 

2 . 3 Linear Modeling 
2.3.1 Equations of Motion 

The controller design procedure begins with the expression of the 
equation of motion in linear time invariant state space form. The non- 
linear, multivariable system that represents the submarine is described 
by: 




f (x ( t ) , u ( t ) ) 



( 2 - 2 ) 



y(t) = g (x(t) ) 



(2-3) 



where 

x is the state vector 

u is the control input vector 

y is the output vector 

These nonlinear equations can be linearized through a Taylor series ex- 

pansion in the vicinity of a nominal point (ideally where — (x(t) ) = 0) 

for small deviations of u(t) and x(t) from the nominal state u , x . 

— — — o — o 

From Eq. (2-2) and (2-3) the linearized dynamics are derived as: 
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Ax(t) s AAx(t) + BAu(t) 
Ay (t) = CAx(t) 



where 



and 



Ax (t) 


= x(t) - 


X 






— o 


Au(t) 


= x(t) " 


u 

— o 




• 

X 






o 


x=x 
o 






u=u 



f (x (t) , u (t) ) 
— o — o 



o 



^(t) 



gCx^t) ) 



3f(x, u) 



A = 



3x 



x=x 
o 



u=u 



( 2 - 4 ) 

( 2 - 5 ) 



( 2 - 6 ) 



( 2 - 7 ) 



B 



3f(x, u) 



3u 



x=x 
o 

u=u 
o 



( 2 - 8 ) 



c 



9g(x) 

3x 



x=x 
o 



( 2 - 9 ) 
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The parametric linearization of the nonlinear equations was accom- 
plished analytically in the interest of accuracy for obtaining good linear 
models, given some of the modeling drawbacks already mentioned. Details 
are given in Appendix A. 

2.3.2 Nominal Point Selection 

A most important step in the design process is to select a linear 
model that best describes the dynamics of the submarine over the widest 
possible range of operating condtions. The control system to be designed 
in this study will use as command variables and, therefore, control, in 
some manner, depth and heading rates. The nominal point chosen for the 
design, then, will be taken about a submarine in a level turn (i.e. , 
pitch angle, or 0 is zero) , since this attitude of a submarine is the 
most likely operating condition. 

All four actuator variables; RPS , <Jr, fib, and <5s must be selected 
to define the particular nominal operating point. The shaft RPS will 
remain at that required for 30 knots in straight ahead motion. The 
rudder deflection, <5r, can be set at arbitrary angles to cause the sub- 
marine to turn at different rates. A rudder deflection of (+) 2 degrees 
was selected for this nominal point design.* 

The next condition involves the use of the stern and bow planes 
that result in zero pitch while turning. For reasons outlined in the 
next subsection, the sail planes will be "locked” at a zero deflection 
angle for this model. With this constraint, the stern planes were per- 
turbed on the nonlinear conputer model in order to achieve zero pitch 
angle. The four actuator variables chosen result in the nominal design 
point of Appendix B. 



The linearization about (-) 2 degrees of rudder would yield different 
system matrices than that obtained about (+) 2 degrees of rudder; this 
is primarily attributed to dynamic asymmetry caused by rotation of the 
vehicle propellor (torque reaction and flow) . 
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2.3.3 Output and Control Variable Formulations 



In this subsection one will gain insight into the appropriate 
selection of the output and control variables to accomplish the proposed 
control design. This will eventually lead to the development, in paral- 
lel, of two controller designs based on two sets of output variables. 

There will be some freedom in the choice of control (actuator) variables 
as well, but their selection will be based on a rather casual approach. 

a. Control Methodology Constraints 

The Loop Transfer Recovery method for the class of Model Based 
Compensators, apart from its advantages, places a very important require- 
ment on the design freedom at a very early stage. It will be seen in 
later sections that the mathematics (singular values et al.) require a 
square system, that is, the number of control inputs be equal to the number 
of output controlled variables. Stated more precisely: 

y (t) = eR m 

u (t) = eR^ 

where R indicates the dimension space of the system. So the requirement 
is 



p = m (2-10) 

and with three independent control surfaces available 

p £ 3 

Note that RPS is not a control variable in the present model. 
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b. Control Input Selection 



It was proposed to have two "y" outputs for the operator to con- 
trol through a " joystick-like " device. From the requirement for a square 
system p = m = 2. How, then, are the two independent control inputs 
selected? Clearly one needs a rudder to produce motion in the horizontal 
plane. For motion in the vertical plane, however, there are two control 
surfaces available, the stern and sail planes, which provide a redundant 
control capability. The stern planes provide about three and a half times 
more force to affect q than the sail planes for the same deflection angle. 
This can be seen by comparing entries B(5, 1) and B(5, 2) of the B_ matrix 
in Table 2.2. Realizing the importance of controlling the stern planes, 
there are two options concerning the use of the sail planes: 

(1) Model the submersible with the sail planes locked at some 

angle, 6b , or; 
o 

(2) Model the submersible so that a linear relationship exists 
between deflection of the stern and sail planes. 

In either method, there is effectively only one independent control sur- 
face for vertical motion control. 

Table 2-2. Model based B matrix. 



:ate 


6s 


6b 


6r 


u 


0 . 29053E-01 


-0 . 15339E-02 


-0 .10980E+00 


V 


0.00000E+00 


0.00000E+00 


0.21864E+01 


w 


-0 . 1860 4E +01 


-0 . 13578E+01 


-0.27856E-05 


P 


0 .00000E+00 


0 . 00000E+00 


0 . 39980E-01 


q 


-0 .45053E-01 


0 . 13042E-01 


0.91875E-06 


T 


0.00000E+00 


0.00000E+00 


-0 . 55462E-01 


$ 


0.00000E+00 


0.00000E+00 


0.00000E+00 


6 


0.00000E+00 


0.00000E+00 


O.OOOOOE+OO 
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For Case 2 above, the sail planes would be scheduled to be de- 



flected as: 



6b 



(-) kSs 



( 2 - 11 ) 



where k is a constant of proportionality. If k is positive, the 
deflection of the sail planes will be in an opposite sense to that of the 
stern planes. This would tend to cause both control surfaces to act on 
q, the pitch moment, in the same direction. For non-zero k, one can 
deduce that less deflection of the stern planes, 6s, would be required 
to produce a given moment q. Since there is a limit on the plane cycle 
rates of about 6 degrees per second, it would appear that have a non-zero 
k could reduce the time lag over a system with constant 6b for a required 

q- 

The revised B matrix, B* , in going from three to two independent 
control inputs according to (2-11) is 



B*(i, 1) = B ( i , 1) - k B (i , 2) 



B*(i, 2) 



B (i , 3) 



( 2 - 12 ) 



where i = 1 to 8 



B(i, 1) = stern plane control vector 



B(i, 2) = sail plane control vector 



and 



B (i , 3) 



rudder control vector 
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In the first option above, a constant deflection other than 6b = 0 
may not prove to be ideal if the controller is to operate at arbritrary 
pitch conditions in a turn. Additionally, if several nominal points were 
selected for gain scheduling, then, when the controller shifted from one 
set of gains to another, a step change in 6b may be required. This change 
could lead to other problems and possible instabilities. 

In this design, k has consequently been chosen to equal zero so 
that 6b = 0. This is consistent with current fleet submarine practice of 
generally not using the sail planes to accomplish high-speed maneuvers. 

The effect of k on robustness and performance will not be addressed in 
this thesis. 

c . Output Variable Selection 

Ideally one would have the present systems directly control the 
bearing rate, i|>, and depth rate z. The expressions: 

\jj = (r cos $ + q sin e)/cos e (2-13) 



and 



z = -u sin 0 + v cos 0 sin $ + w cos 0 cos $ (2-14) 

are nonlinear functions of the state variables. 

To proceed further, the open- loop time domain linear model of the 
plant in Figure 2-2 is introduced. Here the output variable, y (t) , is 
defined as y (t) = Cx(t) . The C matrix will be an important parameter in 
the design of the Model Based Compensator. The C matrix is a constant 
matrix which represents the linearization of (2-13) and (2-14) using (2-9) . 
Letting y^ = jp and y^ = z, the coefficients of the C matrix are given in 
Table 2-3. The matrix indices (1, 2, ..., 8) represent states u, v, w, p, 
q, r, $, 0, respectively, the subscript "o'* indicating the (value of the) 
state at the nominal point. 
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Figure 2-2. Linear open-loop plant model. 



Table 2-3. Linearization of outputs ip and z to obtain the matrix, 
sin (j) 



C(l, 5) = 



cos 0 



cos 



C(l, 6) = 



cos 0 



C(l, 7) = 



-r sin d> + q cos 0 
o o o o 

cos 0 



C(l, 8) = 



(r cos <b + q cos <b ) + tan 0 
o _o o o o 

cos 0 



C(2, 1) 


= - sin 0 

o 




C (2, 2) 


= cos 0 sin 






o 


o 


C (2 , 3) 


= cos <b cos 


0 




o 


o 


C (2 , 7) 


= v cos 0 


cos <(> - w cos 0 sin d> 




o o 


o O O i 


C (2 , 8) 


= -u cos 0 
o o 


- v sin 0 sin (p - w 
o o o o 



26 



To determine how accurately y(t) = Cx(t) represents the actual ob- 

• • 

served outputs ip and z, the author has taken all states of a maneuver of 

a submarine in a turn and reconstructed the linear and nonlinear (actual) 

• « 

representations of \p (t) and z(t). These results are presented in Fig- 
ures 2-3 and 2-4, It is seen that z_ . becomes 50% of z _ . It 

linear actual 

is found that the relative deviation between the linear and actual z 

• • • 

decreases as z . _ increases in magnitude, ip , and ip . are 

actual 3 actual linear 

nearly identical over the full maneuver. The control design based on 

these linearized states will be referred to in the remainder of this 

• • 

study as the ip - z controller. 

An alternative to a C matrix which is dependent on the vehicle 
states at the nominal point is to choose such that its entries are 
independent constants , but that would still represent rates that are 
desirable to control. 

Referring to Figure 2-4 once more one can see that such a term 
for ip exists. The state r has , throughout this maneuver , remained within 
4% of ^ actua ]_- This is made clear by referring to (2-13) and noting 
that for nearly all turning maneuvers, r >> q and the cos 0 = 1. To 
control depth rate, one can revert back to current submersible opera- 
tional practice of commanding pitch angle. This follows from (2-14) 
where u >> v, w so that z = -u • sin 0 = -u0 . Since forward velocity is 
known, z is proportional to 0 so C becomes: 

C(l, 6) = 1 

C(l, 8) = 1 

The designed based on this C matrix will be known as the r - 0 controller. 

• • 

For implementation purposes the exact values of ip and z are always 
available (observable) for both the nonlinear models and the actual submer- 
sible. Observer devices to measure z include depth gauge measurement. For 
more accurate designs with less environmental noise, gimballed accelerom- 
eters could be integrated to provide z. The installed rate gyros used 
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- (ft/s) 




Figure 2-3. Reconstruction of actual and linearized depth rate for 
the submarine in a turn. 
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deg/s 




Figure 2.4. Reconstruction of actual and linearized turn rate, and the 
state "r" for the submarine in a turn. 
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for fire control system can measure the submarined true heading rate, ip , 
regardless of the vehicle attitude. The r and 0 states can be measured 
with similar ease. 

2.4 Dynamic Analysis 

In this subsection the open-loop poles (eigenvalues) and eigen- 
vectors will be calculated, and the conditions for controllability and 
observability of linear control systems will be presented. 

The eigenvalues of the A matrix resulting from the linearization 
are shown in Table 2-4. The corresponding normalized eigenvectors are 
shown in Figure 2-5. It would appear that all modes are dominated by 
the velocity states u, v, and w. A different picture results when the 
values of the eigenvectors involving angle and angle rate states are 
scaled from radians to degrees and the velocities u, v, and w remain 
in feet/second. Although the choice of units was somewhat arbitrary, 
the scaled and normalized eigenvectors of Figure 2-6 seem to provide a 
clearer picture of the state contribution to the system modes. The com- 
plex pole, -0.19 ± 0 . 32i , is associated with the natural roll behavior of 
the submarine. The slowest pole, -0.013, is primarily a pitch mode. The 
zeroes of the system will be discussed in Chapter 4, after the addition of 
augmented dynamics. 

Table 2-4. Complex eigenvalues of A. 



-0 . 190874E+00 
-0 . 190874E+00 
-0 . 487435E+00 
-0.443184E+00 
-0 . 978560E-01 
-0 . 432890E-01 
-0 . 591569E-01 
-0 . 133116E-01 



0.320975 I 
-0.320975 I 
0.000000 I 
0.000000 I 
0.000000 I 
0.000000 I 
0.000000 I 
0.000000 I 
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Figure 2-5. Unsealed eigenvectors of [A], 
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Figure 2-6. Sclaed eigenvectors of [A] . 
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Controllability and observability properties need to be established 
to proceed further with the control system design. Controlability is en- 
sured only if the matrix 



[B, AB, 




(2-15) 



has rank n. Similarly, observability is guaranteed if the matrix 



CA 




(2-16) 



has rank n. 

A simpler, but more conservative method than (2-15) and (2-16) 
above will be used to establish controllability and observability. Since 
eigenvalues of A are distinct and nonzero, one can find the complex modal 
matrix, T, whose columns are the distinct eigenvectors, v_^. For any set 
of scalars c. (not all zero): 

l 



2 C.v = 0 

i 1 

implies that T is a nonsingular matrix. Then the state equation can be 
written in the modal domain [5] : 

d -1 -1 

— x* = T ATx* + T Bu 
dt — — — 

y = C T x* 
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The controllability and observability conditions are then more than 
satisfied by observing that the real parts of the columns of T and 

the rows CT have no zeros. This was found to be the case for both the 

• • 

r - 0 and z - ip systems. 

2 . 5 Performance Specifications 

Performance specifications outlined in this section are not all 
encompassing and may not coincide with established Navy specifications 
for submersible control. The performance requirements as stated here 
are mainly driven by the intuitive engineering approach to obtain good 
command following, speed of response, robustness, and disturbance re- 
jection with due respect to the natural dynamics of the vehicle as 
they emerge from the available model. These performance requirements 
(guidelines) will be accomplished through loopshaping techniques to be 
discussed in Chapter 3. 

Since there are no natural integrator states (the 8 x 8 A matrix 
is nonsingular) , elimination of steady-state error to step inputs is not 
possible. Good command following and elimination of steady-state errors 
will be obtained through the use of integrators in the command variable 
channel. Although an actual rate control system may also experience 
ramp-like commands from the "joy-stick", the command inputs here are 
closer to being a step due to the large time constants of the system. 
Therefore, the system described in this thesis will not be designed to 
meet type-2 specifications. 

The open- loop dynamic simulation demonstrates that achievable 
vehicle settling times due to deflections of control surfaces are a 
function of the vehicle speed. For a submersible at a speed of 30 knots, 
settling time of 50 - 60 seconds are achievable. When settling is de- 
fined to occur in 4 or 5 time constants, the resulting minimum band- 
width requirements for crossover are 0.07 to 0.1 radians/second. From 
the performance aspect it is desirable to have all channels crossover 
at roughly the same frequency. 
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On the high frequency side, one must be able to reject noise and 
possible modeling errors. Noise sources generally originate from the 
environment (true sensor readings) or from the sensor itself. Sensor 
noise typically comes at a higher frequency than the system bandwidth and 
should not affect ship dynamics since ship eigenvalues will typically lie 
in the lower frequency band. It is desired that the ship actuators do 
not respond to normal environmental disturbances that may cause unnces- 
sary actuator motion. Of concern here is the effect of the submersible 
operating close to the surface, but not so close as to be in the danger 
of broaching, and experiencing excitations of surface waves. Based on 
Figure 2-7 [6] , it is apparent that a typical wave spectrum has a fre- 
quency range between 0.2 and 2 radians/second. The excitation or driving 
frequency the submarine experiences will be 



where 

03 

e 

0) 

k 

u 

$ 



a) = a) + ku cos ($) 
e 

frequency of encounter 

wave spectrum frequency 

oo^/g = wave number in deep water 

ship forward velocity 

ship head relative to sea direction 



then, based on an estimated max ship speed of 35 knots (60 fps) : 



0 ) max = 9.5 rad/s (in head seas) 

e 



0 ) min 
e 



0.13 rad/s (in following seas) 



35 




(rad/sec) 



Figure 2-7. Typical form of a wave spectrum 

containing swell (from Ref. [6]). 

and therefore at any vehicle speed, 35 knots or less, all other ship 

headings relative to the sea will be contained within oj max > co > co min. 

e e e 

The proposed controller will, then, be able to meet these minimum 
criteria with the suggested bandwidth, at least in the region where linear 
behavior dominates the vehicle dynamics. 
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CHAPTER 3 



CONTROL SYSTEM DESIGN METHODOLOGY 



3.1 Robustness 



3.1.1 Robustness Design Considerations 

Control system design frequently involves the tradeoff of realizing 
high performance while preserving good command following, disturbance 
rejection, and system stability. The achievement of high performance is 
pointless if the control system becomes unstable as a result. This is 
particularly important in the present design. Overall system instability 
could cause the submersible, in some instances, to exceed its safe operating 
depth or ground itself, resulting in the loss of the submarine and its crew. 
The multivariable methods used to assure stability and robustness will be 
reviewed next. 



system stability and robustness in the presence of unstructured uncertainty. 
The singular values of an (n x n) complex matrix A derived from the spectral 



Singular value analysis has evolved as a reliable method to evaluate 



norm 



| | a[ j 2 are defined by: 




i = 1, 2, 



. . . , n 



l 



l 



where 



A is the complex conjugate transpose of A 



is the ith eigenvalue operator 
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and represents a measure of how close the matrix A is to being singular. 

A plot of the singular values of a multivariable transfer matrix could be 
interpreted in an analogous manner to the single input — single output Bode 
plots. 

The properties of singular values will now be applied to obtain 
statement of guaranteed stability. Let Figure 3-1 represent the MIMO 
feedback conf iguration , where G(s) is the true plant transfer function at 
all points. Figure 3-2 represents the same system except now G(s) = 
G(s)L(s) , so that L(s) represents a perturbation of the nominal system 
G(s) . It is desirable, then, to determine under what conditions of L(s) 
the nominal system will remain stable. 



£-Q 

a 



Figure 3-1. MIMO feedback conf iguration. 



K(s) 



G(s) 



y(t) 




Figure 3-2. MIMO feedback configuration representing 
perturbation of the nominal system. 

There are many ways to represent L(s) . To characterize the robust- 
ness of the system, Lehtomaki [7] utilizes matrix transfer functions located 
at various points in the plant. He defines addition, multiplication, 
division, and subtraction errors. Multiplicative error has the property 
of preserving relative quantities and is independent of the units used 
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in modeling the system dynamics and therefore will be used in this model. 
Furthermore, a premultiplicative error of the loop transfer function will 
be defined to reflect errors at the plant output, as shown in Figure 3-3. 
In this instance 



L (s) = I + E (s) 



so that the premultiplicative error, E (s) , is found to be 

— pre 



E (s) 
—pre 



= [G (s) - G(s)]G _1 (s) 



(3-1) 




Figure 3-3. MIMO feedback configuration reflecting 

premultiplicative error at plant output. 



From the Fundamental Robustness Theorem and the Theorem for Robust- 
ness for Multiplicative Errors (see Lehtomaki [7]), for guaranteed stabil- 
ity, the following inequality must hold: 



a [E ( jw) ] < a . [I + (G(jw)K(ja)) ) 1 ] 
max —pre J — mm — — — 



(3-2) 
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If the inequality (3-2) above is not true, it does not necessarily 
imply the system must be unstable. Singular value analysis does not take 
into account the actual direction of the perturbation matrix that may 
cause a system to become singular; hence the system may actually be more 
robust than singular value analysis may indicate. Therefore the inequality 
(3-2) above is conservative. 

3.1.2 Scaling 

Scaling is a method of weighting the physical units of a system, 
through an appropriate transformation, so that the numerical values of 
the variables become equally significant. Scaling and its effect in 
designing robust multivariable control systems has been resently dis- 
cussed by Kappos [9] and Boettcher [8] . Apparently, scaling does in 
fact change the singular value magnitudes, but it may not necessarily 
change the robustness of the system [8]. Ostensibly, not all effects of 
scaling are understood at present. 

A systematic method does not yet exist to obtain optimal scaling. 
Kappos [9] selects a scaling matrix based on the expected nominal departure 
(error) of the output. Kwakernaack [10] , although not specifically addres- 
sing the robustness issue, similarly suggests weighting the states by 
their tolerable error, i.e. , a deviation of 10 ft/s in velocity may be 
as bad as a deviation of 0.2 radians in pitch, and then base the weighting 
on the ratio of the tolerable errors. In this thesis the author presents 
a method that allows the system to scale itself through loop shaping 
techniques. This procedure will be presented in Chapter 4. 

A scaling transformation to transform state vectors must be diagonal 
and positive to preserve the system eigenvalues. Briefly, if u, y, and x 
are the original state vectors, and u* , y 1 , and x* are a set of scaled 
vectors, then: 
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u = S u' 



Y = S y ' 



x = S x' 
— — x— 



(3-3) 



consequently for a system of the original form: 



x = Ax + Bu 



y = Cx 



(3-4) 



the transformed matrices become 



A' = S 1 A S 
— —x x 



B' = S 1 C S 
— — x y 



C' = s ' 1 c s 

- -y X 



(3-5) 



Scaling directly affects transfer functions as well. Defining the 
open-loop transfer function, G^fs) °f Figure 2-2 as 



Gq L (s) = C(sl - A)" 1 B 



(3-6) 



it can be shown straightforwardly that the scaled transfer function be- 
comes 



G' (s) = S -1 G (s) S 

— — y — — i 



(3-7) 
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Similarly, scaling E (s) , the premultiplicative error defined in Eq. (3-1) , 

pre 

one obtains 



E'(s) = S _1 E (s) S (3-8) 

_ -y _ -y 

Finally, a useful property of the scaling transformation will be 
developed next. Take the special case of Eq. (3-7) where is the 

identity matrix so that 



G’ (s) = 



S 1 G(s) 

-y ~ 



(3-9) 



By taking the complex conjugate transpose of Eq. (3-9) above: 



G' H <s> - G H (s) S’ 1H 

— — -y 



(3-10) 



From Eq. (3-9) and (3-10) then 



G'(s)G ,H (s) = 



— 1 H — 1H 

S G (s) G (s) S 

-y - - -y 



(3-11) 



and, applying the two- norm implies 



G' (s) 



l ^ 1 Skills 



(3-12) 



so that 



G' (s) 



Is ; 1 £(s) | i 2 



(3-13) 
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Note that S is not a function of s. Through matrix multiplication, 

-1 -1 -1 

the diagonal entries of S 1 . .. S of the matrix S independently multiply 
the rows of G(s) . Knowing the following relation to hold 

|aA| | = | a | • | |a| 

where a is a scalar , 
it follows that the i singular values of | [s^ 1 G(s) | | satisfy 

CT.[G'(s)] = sT 1 a . [G (s) ] (3-14) 

1 — 11 — 

Hence, it is seen that the transformation matrix S ^ directly effects the 

Y 

magnitude of the singular values of the original transfer matrix. Equation 
(3-14) will be an important tool that will be employed in Chapter 4 to help 
select appropriate scaling transformations to achieve consistent loop 
shaping. 

3 . 2 LQG Compensator 

The Model Based Compensator (MBC) has evolved from the optimal esti- 
mation theory (Kalman Filter) and optimal control theory (Linear -Quadratic 
Feedback) . The concepts presented here can be found in more complete form 
in references [10] and [11] . The form of the compensator used in this 
design together with the state definitions is shown in Figure 3-4 with 
the plant and compensator dynamics represented separately. The transfer 
function definitions used in the remainder of this paper will be based on 
Figure 3-4. 

Referring to Figure 3-4, the error e(t) is the input and u(t) is 
the output of the compensator; the overall MBC transfer matrix is defined 
by 



u(s) = K ( s) e (s) 
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Figure 3.4, The model based compensator in a feedback configuration. 



where 



K(s) = G (sjl - A + BG + HC) 1 H (3-15) 

Also, it is easily shown that the overall closed- loop transfer function, 
(s) is 



2cl < s) 



det(sl - A + BG) det(sl - A + HC ) 



(3-16) 



It follows that the requirements for stability of the system are: 



Re " BGl < 0 



(3-17) 



and 



R e A i [A - HC] < 0 (3-18) 

which are the poles of the compensator. 

The conditions for ensuring a stable compensator are now postulated. 
From linear system theory, if the pair [A, B] is controllable (or stabiliz- 
able) , then there exists at least one feedback gain matrix G to ensure that 
all closed-loop poles of [A - BG ] are in the left half s-plane. G is 
given by 



G 



-1 



R 



BK 



(3-19) 



where K solves the Control Algebraic Ricatti Equation (CARE) : 



T -IT 

0 = - KA - A K - Q + KBR B K 



(3-20) 
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Likewise, if the pair [A, C] is controllable (or detectable), then 
there exists at least one output gain matrix H_ such that all closed- loop 
poles are in the left half s-plane. The matrix E is given by: 



H 



„ T -1 
E C = 



(3-21) 



where E_ solves the Filter Algebraic Ricatti Equation (FARE) : 

0 = EA + A T Z_ + LHL - ZCV 1 ^ (3-22) 

3 . 3 Loop-Transfer Recovery 

In the sequel, some properties of the Kalman Filter that forms 

the basis for the Loop-Transfer Recovery (LTR) method to be employed in 

this design will be discussed. The Kalman Filter loop transfer function, 

G. is taken at point 1 on Figure 3-4. The resulting loop-transfer func- 
KF 

tion is defined 



G = C(sl - A) 1 H (3-23) 

— KF — — — — 

One can also define G^^, the filter open- loop transfer function as 

G„ t = C (si - A) -1 L (3-24) 

rOL 

where L is a free design parameter that is chosen to give G desirable 

— — rCJL 

singular values to meet the performance and robustness specifications. 
The Kalman Filter Domain Equality 

[I ♦ 0^(3) HI + 6^.(8) ] H - I G F01 (S) G» 0 l< S ) (3-25) 

has been derived from Eq. (3-22) , (3-23) , and (3-24) where the substitu- 
tion 
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_0 = yl for ]i > 0 



has been made in the FARE Eq. (3-22) . 



As 



in Eq. (3-25) , it follows directly from the definition of singular values 
that 



which is valid for low frequency range, and represents the recovery of the 

loop shapes of the G from G . The value of y is also a free-design 

KF FOL 

parame ter . 

The Kwakernaack Loop Transfer Recovery process (see e.g. , Doyle, 
Stein [1], Kwakernaack [12]) is described next, in which the shapes of 
the Kalman Filter singular values of Eq. (3-26) are recovered in the loop- 
transfer matrix of the Model Based Compensator. The loop- transfer matrix 
of Figure 3-4 is given by 




(3-26) 



T = G ( s) K(s) 



C (si -A) 1 BG(sI - A + BG + HC) 1 H 



(3-27) 



The LTR method sets 



R 



I 



and 



T 

Q = qC C 



47 



in the CARE Eq. (3-20) . Then, as 



the corresponding in Eq. (3-27) becomes 'large' so that Eq. (3-27) 
becomes 



for low frequencies. 

The LTR method is guaranteed to work in the open- loop system is 
minimum phase. For non-minimum phase systems there is no such guarantee, 
although some recovery of performance and robustness properties are 
expected as the non-minimum phase zeros move further away from the desired 
operating bandwidth (see [13]). 

Robustness Theorems 

The following robustness properties can be derived from the Kalman 
Filter Equality (3-25) and are required for system robustness: 



T ro = C(sl - A) -1 H 



(3-28) 



and, consequently 




(3-29) 




(3-30) 



which follows directly from the KFE ; and 



( 2 ) 




(3-31) 



which has been proved by Kappos [9] . 
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CHAPTER 4 



COMPENSATOR DESIGN PROCEDURE 



4. 1 Augmented Dynamics 

Augmenting the dynamics of the submersible control system serves a 
dual purpose. One is to model the actuator dynamics to make the model as 
accurate as possible, and achieve desirable roll-off at crossover for 
robustness. The other is to include integrators to cause the compensator 
to behave as a type-1 system, which will permit the submersible to achieve 
zero steady- state error to step inputs and disturbances (i.e. , good command 
following) . The result will be to increase the order of the sytem by four 
in the present two- input design. 

A block diagram of the augmented model appears in Figure 4-1. It 

is seen that the augmented dynamics have been placed in the command channel. 

u is the true commanded input (a physical variable), u , the output of 
c act 

the augmented dynamics, is not. The mathematics of the augmented states 
will be manipulated in such a way as to provide a means to achieve the 
desired loop shapes of G . 



Actuator Dynamics 

The complete actuator dynamics are governed by their mass proper- 
ties of the rudder and planes, and angle rate limits imposed by the 
electro- hydraulic servomechanisms that position them. Since the rate 
limits are highly nonlinear, they will not be modeled per se . Instead, 
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AUGMENTED DYNAMICS EXISTING PLANT 

Figure 4-1. Augmented dynamics placed in command variable 
channel of plant. 

the compensator will be designed so that the rudder and stern planes are 
never driven past their rate limits of 7 and 6 degrees/second respectively, 
nor past their maximum allowable deflection angle. 

The actuator dynamics of the rudder are represented by a second- 
order system. These are given for the stern planes as 

oj- = 1.37 

d 

C = 0.9 



and for the rudder as 



ok = 1.09 

d 

e = 0.9 

as C approaches unity, a second-order system begins to lose its oscillatory 
characteristics and behavior tends toward a first-order system. Since this 
is very nearly an over damped system, the actuators will be modeled as a 
first-order system. By Ogata [15] the rise time for a second-order system 
is given by 
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t 

r 



tt - tan 



-l A - c‘ 



so that for the stern planes 



and for the rudder 



t = 1.96 s 

r 



t = 2 . 47 s 

r 



Equating the rise time to 4 time constants in an equivalent first-order 
system , then 



'6s 



= 0.49 s 



and 



"6r 



= 0.62s 



The augmented actuator dynamics, with integrators, are presented in 
Figures 4-2 and 4-3. The full state space representation of the dynamics are 



I- 1 



{x a } - 



6s 

0 



0 



' 6r 



6s 

0 



0 



* 6r 

0 

0 



(x } + 

a 



o o 



0 0 



1 0 



{u} 



and 



10 0 0 



0 10 0 



{x } 
—a 



C x 
—a— a 
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Figure 4-2. Stern plane command channel with 
augmented dynamics. 




Figure 4-3. Rudder command channel with 
augmented dynamics. 

These dynamics are introduced to the 8th order state space repre- 
sentation to produce a 12th order system. To get it in a form for further 
design manipulation a new state vector is defined 




where x* is the unaugmented 8th order system augmented with the two command 

channel actuator time lag states and x are the two integrator states in 

— s 

the command channel. The augmented system can be written as 
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Y 




Cx 



(4-2) 



This form of state space equations will be used in Section 4.2 to derive 
suitable loop-shaping algorithms. The A and 13 matrices of the augmented 
system are found in Appendix B. 



Augmented System Dynamics 

The four additional states add two poles at the origin and one each 

at -1.6 and -2.0 to those already listed in Table 2-4. The multivariable 

• • 

zeros for both the r - 9 and ^ - z system are listed in Table 4-1. Note 
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1C° 




Figure 4-4. Unaugmented open- loop plant; r 
singular values vs frequency. 




0 system; 
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Sl'jrjlll/IR VALUES (LOG MHGN 1 TUDE ) 

• 2.50 1.25 0.00 1.25 2.50 3.75 



OPEN LOOP PLRNT 2X8 0Z 




Figure 4-5. 



Unaugmented open- loop plant; if) 
singular values vs frequency. 



z 



system; 
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SlNGULflR VALUE'S (LOG MAGNITUDE) 

10.0 7.5 - 5.0 - 2.5 0.0 2.5 



OPEN LOOP PLRNT 2X12 R, THETR UN5CRLED 




\ 



n i ! 1 ' — r- 



10‘ 1 10 
FREQUENCY (PP0/5EC) 



10 



Figure 4-6. Augmented open- loop plant; r - 0 system; 
singular values vs frequency. 
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Table 4-1. Augmented system zeros (finite) . 



r - 9 System Zeros: 



NUMBER 



DAMPING 



FREQUENCY 



1 -5.796252290E-02 

2 -1.570545617E-01 

5 -2. 9Q7626130E-01 

5 -1. 8820000** 2E— 01 -5 .060423035E-01 5. 2S8274395E-01 



3.592786288E-01 



^ - z System Zeros : 



NUM3ER REAL IhAG DAMPING 

1 -5.796476182E-02 

2 -2 . 38 779884 7E-01 

5 -5.115217184E-01 

5 -1.813410260E— 01 -5. 0882406 98E -01 S.063557471E-01 

6 7.209455517E-01 



FREQUENCY 



5.581296886E-01 



• • 

that ^ — z system has a non-minimum phase zero at 0.7, the existence of 
which is most likely due to the relationship between the pitch (0) and 
heave (w) states in z. It is pointed out here that both the actuator 
dynamics and the non-minimum phase zero are beyond the desired system 
bandwidth of 0.1. Therefore, they are not expected to effect the design 
appreciably. However, the actuators have been incorporated in the design 
to provide additional roll-off at higher frequencies. 

Open-Loop Transfer Function 

The singular values of the unaugmented (8th order) G systems are 

OL 

plotted in Figure 4-4 and 4-5. In both of these figures, the minimum 
singular values correspond primarily to the turn related outputs, and r 
The beginning of gain roll-off corresponds to a system pole as it would 

in a Bode plot. The maximum singular values in both figures begin to 

-2 

roll-off at 10 radians/s, which corresponds to the pole at -0.013 
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51UGULflR VALUES (LOG MAGNITUDE ) 

• 5.0 2.5 0.0 2.5 



OPEN LOOP PLANT 2X12 ZDOT, PSIDOT UNSCALED 




• • 

Figure 4-7. Augmented open- loop plant; \j) - z system; 
singular values vs frequency. 
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(the 0 - w dominated eigenvector of Figure 2-6) . The minimum singular 
values start to roll-off at 10 ^ radians/s, which relates to the pole at 
-0.097. This pole corresponds to the "slowest” eigenvector to first con- 
tain components of the r state. The bandwidth corresponding to the r and 
t]/ outputs is nearly an order of magnitude greater than the z and 0 outputs. 

The open-loop singular values of the augmented (12th order) system 
are shown in Figure 4-6 and 4-7. The integrators in the command channel 
correspond to the slope of 1/s or 20 dB/decade at low frequency. Also note 
that the minimum singular values of both figures are nearly identical over 
the chosen frequency range. It will be shown that either choice of r or ^ 
as an output variable will result in nearly equal performance in the con- 
troller design. 



4.2 Loop-Shaping Techniques 

In Chapter 3 the notion of the free-design parameter p and matrix L 
was introduced to provide a means of shaping the singular values to meet 
the design specifications. In this section, the author will derive a pro- 
cedure to obtain a. [G ^ ] that are tied together throughout the desired 
l — FOL 

operating bandwidth, and will eventually result in a desirable tight cross 
pattern of the Model Based Compensator loop transfer function T (s) . 

The filter open- loop transfer function is given in Eq. (3-24) as 



« fol - £<=!-« £ 



Let L be of the form 






L = 



—2 
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where 



is all states excluding the integrator states 
is the integrator states 



and, starting with the full augmented system dynamics in the state space 
form as developed in Eq. (4-2) 



A 





one obtains 



si - A 




-B 



si 



Taking the block inverse 



sl_ - A) 



-1 



( S I - A 10 } 



-1 



0 



(si 





s 




(4-3) 



Now for low frequencies 



s -> 0 
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-a’ 1 

-10 



-1 

-A B 
“10 — I 



h. 



C (Si - A)” 1 L 



[ ^10 



0] 



S 



l_r*2 



, -c a -H. £ ^ 
—10— 10^-1 s 



Hence, for small s, Eq. (4-4) further becomes 



(4-4) 



C (si - A)" 1 L s 



- £ iio ~t h 



(4-5) 



By Eq. (3-22) 



— FOL ( S } 



= C (si - A) 



-1 



To tie the singular values together, it is desirable to set 

I 



-FOL a s 



and, consequently, obtain from (4-5) 



- £ 4o 2, 



(4-6) 



Finally, solving for L 2 



±2 ' - 



(4-7) 
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The dependence of the low frequency singular values should be clear; 

]±2 drives the controller through the integrator states. 

An analogous method was used to the high frequency singular values 
together. From Eq. (4-3) it follows that 

T T -1 

k ± = C (CC ) (4-8) 

4. 3 Model Based Compensator Design 

Equations (4-7) and (4-8) were applied to both designs. Figures 

4-8 and 4-9 are the plots of the shaped a. (G„ ) which show how the 

l FOL 

singular values are tied together at low and high frequencies. By varying 
the design parameter y , these singular values can be shifted up or down. 
However , neither plot would result in a tight crossover pattern at the 
desired crossover frequency of 0.08 radians/s. It is through the use of 
the scaling transformation properties developed in Eq. (3-14) that the 
shapes of these singular values approach the desired shapes for good com- 
mand following, disturbance rejection, and crossover frequency specifications. 

Referring to Figure 4-8 it is proposed that the minimum singular 
value, r, be scaled to match the maximum singular value, 9, at its greatest 
separation. This point occurs at a separation of about 24 dB, which is 
equivalent to a ratio between the values of r and 9 of 16, at a frequency 
of 0.02 radians/s. For numerical convenience, r will be scaled by a factor 
of 10 (versus 16). This means that r will be observed in "tens" of radians, 
while the output 9 continues to have units of radians. 

A similar approach is applied to figure 4-9. In this case, the maxi- 
mum separation approaches 4 orders of magnitude. For numerical convenience, 

\ p will be sclaed up by a factor of 10 (the same factor used for r above) , 
where z is scaled down by 0.01. 
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SINGULAR VALUES (LOG MAGNITUDE) 

- 2.50 - 1.25 0.00 1.25 2.50 3.75 



PLANT 2X12 CtSI-flf'l R,THETR, UNSCALED, SHAPED 

C 




Figure 4-8. Shaped and unsealed filter open- loop transfer 
function? r - 0 system? singular values vs 
frequency. 
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SINGULAR VALUES (LOG MAGNITUDE) 

-2.50 M.25 0.00 1.25 2.50 3.75 5.00 



PLANT 2X12 CfSI-Rl 1 L Z,PSI DOT, UNSCALED 




Figure 4-9. Shaped and . unsealed filter open- loop transfer 
function; ip - z system; singular values vs 
frequency. 
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Consequently the transformation matrices, 




become 



10 0 “ 

0 1 _ 

for the r - 0 system and 

10 0 
0 0 . 01 . 

• • 

for the ip - z system. 

Figure 4-10 and 4-11 show the effect of scaling alone on the open- 
loop transfer function C[sl_ - A] 1 B. The transformations of (4-9) and 
(4-10) have resulted in the open-loop transfer function being scaled up or 
down by the scaling factor applied to each output as predicted in Eq. (3-14) . 
Note that in these plots the maximum and minimum singular values are plotted 
instead of the singular values of the outputs. This accounts for the appear- 
ance of the singular values of Figures 4-10 and 4-11 of being pinched to- 
gether. If a. (G ) had been plotted, the singular values would have inter- 
i OL 

sected at this point. 

The final loop shapes of G are found in Figure 4-12 and 4-13 where 

— FOL 

both the scaling transformation (3-9) and loop shaping Eq. (4-7) and (4-8) 
have been employed. These loop shapes are nearly ideal for ensuring a tight 
crossover pattern and nearly 20 dB/decade roll-off throughout the operating 
frequency range. The current crossover frequency of 1 radian/s is to be 
expected due to the asymptotic behavior of a^(I/s). 

The final parameter to be chosen is y. By Eq. (3-26) 






°i[^ W j “’] 
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SIfICULFiR VALUES (LOG MAGNITUDE) 

- 10.0 7.5 5.0 - 2.5 0.0 2.5 5.0 



OPEN LOOP PLRNT 2X12 R, THETfl SCALED 




Figure 4-10. Scaled open- loop plant; r -0 system; 
singular values vs frequency. 
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SINGULAR VALUES (LUG MAGNITUDE) 

- 5.0 - 2.5 0.0 2.5 ! 



OPEN LOOP PLANT 2X12 ZDOT, PSIDOT SCALED 







• • 

Figure 4-11. Scaled open- loop plant; \p - z system; 
singular values vs frequency. 
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SINGULAR VALUES (LOG MAGNITUDE) 

- 2.50 1.25 0.00 1.25 2.50 3.75 5.00 



PLANT 2X12 C(SI-fl) H L R, THETA, SCALED, SHAPED 




Figure 4-12. Shaped and scaled filter open- loop transfer 
function; r - 0 system; singular values vs 
frequency. 
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SINGULAR VALUES (LOG MAGNITUDE) 

- 2.50 - 1.25 0.00 1.25 2.50 3.75 5.00 



PLANT 2X12 C(SI-A) H l Z,PSI DOT, SCALED, SH4PED 




Figure 4-13. Shaped and # scaled filter open- loop transfer 
function; ip - z system; singular values vs 
frequency. 
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and solving for y then 



= {a [G (a) ) ] } 
max — : FOL o 



-2 



where 5^ is the (maximum) desired crossover frequency. The value of 
obtained were 126 and 144 respectively, for the r - 0 and jp - z controller, 
using a maximum crossover frequency of 0.09 radians/s. 

With the design parameters, y, L, and s” 1 thus determined, the design 

- -y 

can proceed with the application of the FARE Eq. (3-22) and (3-21) to 
calculate the Kalman Filter gain matrix H. The singular values of 
are plotted in Figures 4-14 and 4-15. The loop shapes have been re- 
covered in both cases with crossover in the range of 0.07 to 0.09 radians/s. 

Finally, the controller gain matrix, K is calculated by the CARE 
Eq. (3-20) and (3-19) . The only free design parameter to choose is q, the 
control weighting index. Generally, to prevent control saturation from 
high gain of K, it makes sense to choose q as small as possible and still 
meet the design specifications. In both controller designs, when q was set 
equal to 10 , a crossover of 0.075 radians/s was obtained by the singular 
values of the loop-transfer matrix, T(s) . These are plotted in Figures 4-16 
and 4-17. Note the roll-off at crossover is 20 dB/decade , as expected, 
and that Eq. (3-27) fails to hold at frequencies above 0.7 radian/s. 



Robustness 

The singular value plots of the Kalman Filter return difference and 

Kalman Filter inverse return difference matrices are shown in Figure 4-18 

• • 

and 4-19 for the r - 0 system, and in Figure 4-20 and 4-21 for the - z 
system. The robustness criteria set forth in Eq. (3-30) and (3-31) have 
been met. 

A discussion of the robustness of the loop-transfer matrix according 
to Eq. (3-2) is given in Chapter 5. 
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INI? 1 '1. MR VFILUCS (LOG MAGNITUDE) 



PLANT 2X12 C(SI-A)-1H 




1 i — i — r~' i ~f in 1 1 — ; — i — r-r-n 

:o ' 3 iz* io ' 1 :o fl 

FREQUENCY [RflD/SEO 



Figure 4-14. Kalman filter loop-transfer function; r - 0 system; 
singular values vs frequency. 
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SI I JGULFIR VALUES l LOG MflGN I TUDEi ) 

- 2.0 - 1.0 0.0 1.0 2.0 



PLANT 2X12 C(S1-A)-1H, Z,PSI 




1C" 



10 "' 



10 

FREQUENCY 



10 

(RAD/SEC) 



10 



10 



Figure 4-15. 



Kalman filter loop transfer function; ip 
singular values vs frequency. 



z 



system; 
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SINGULAR VDLUEIS (LOG MAGNITUDE:) 

• 12.5 - 10.0 - 7.5 - 5.0 - 2.5 0.0 2.5 




FREQUENCY (RRD/SEC) 

Figure 4-16. Loop-transfer recovery; r - 0 system; singular 
values vs frequency. 
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SINGULAR 1/flLUES (LOG MAGNITUDC) 

10.0 8.0 - 6.0 4.0 2.0 0.0 2.0 



COMPENSATED SYSTEM 30R2 0=10 Z, PSI DOT 




« • 

Figure 4-17, Loop-transfer recovery; ip - z system; singular 
values vs frequency. 
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SINGULAR VALUES (LOG MAGNITUDE) 
0.8 1.2 1.6 2.0 



( 1 + G) R-THETR COMPENSATOR 




Figure 4-18. Kalman filter return difference; r - 0 system; 
singular values vs frequency. 
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Mfi'-Ji r-; values nw; mGMiTitDi:) 



!T + G-n R-THETR COMPENSATOR 
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I 



I 




Figure 4-19. Kalman filter inverse return difference; r - 0 system; 
singular values vs frequency. 
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SINGUU1R vm.UCS (LOG IKlGNITUDd 



a + G) PS I - 2 DOT COKPENSfilOR 




Figure 4-20. Kalman filter return difference; \p - z system; 
singular values vs frequency. 
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Figure 4-21. 



Kalman filter inverse return difference; ijj 
singular values vs frequency. 
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system? 
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CHAPTER 5 



MODEL BASED COMPENSATOR EVALUATION 



5 . 1 Compensator Design Implementation 

In this section, the performance of both the r - 0 and ip - z con- 
trollers will be tested on the nonlinear submarine simulation to see how 
closely the performance specifications are met and to check for instabili- 
ties. Note that the linear plant dynamics G(s) presented in the feedback 
configuration of Figure 3-4 has been replaced by the nonlinear plant 
dynamics of the simulation. The states at point 2 of Figure 3.4 are 
therefore the actual measure of output variables. This means that the 
error vector, £(t) , at the input of the Model Based Compensator is always 
the true difference between the commanded input and the output variables. 

In order to avoid the additional computations required to transform 
the output variables back to their original units, the output variables are 
handled in their scaled form. To maintain a properly scaled error vector, 
e_(t) , the true output variables and the command input variables, £(t) , are 
multiplied by the appropriate scaling transformation. Figure 5-1 illustrates 
how the model based compensator feedback configuration for this design imple- 
mentation is modified by scaling. 

5 . 2 Evaluation of Output Variable Selection 

Comparison of the two system performances will be accomplished in this 
section to determine which rate control system provides better performance 
and robustness characteristics. The performance evaluation will be based 
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Figure 5-1. Modifications of the MBC feedback configured 
due to scaling. 

on commanding various bearing rates, while ordering zero pitch or zero 
depth rate to the two compensators. The resulting behavior of the vehicle 
with respect to depth will then provide a common means of comparing the 
response of the two designs. 

The first simulation run for both controllers is for a small ex- 
cursion about the nominal point of 1 degree/s and zero feet/s (or zero 
pitch for the r - 0 system) ; it is expected that the nonlinear model 
should behave linearly in this vicinity and would provide a gauge of 

how successful the model was in meeting some of the performance specif i- 

• • 

cations. The entire state outputs for the ip - z compensator simulation 
are found in Figure 5-2. The ip and z command errors under the control 
input heading of Figure 5-2 indicate the performance. It can be seen 
that settling times for both outputs are indeed within 50 - 60 seconds 
and are in compliance with the performance specifications of Section 2.6. 
Note also that there is already a small error in z after 60 seconds. 

This error tends to zero as more time is allowed for the simulation, and 
results in a maximum vehicle depth change of about 5 feet for this maneuver. 
The command errors of the r - 0 compensator shown in Figure 5-3 exhibit 
nearly the same performance as the ip - z compensator for the heading rate 
orders, but the 9 output state experiences a small overshoot. Figure 5-3 
also shows that as 0 approaches the value of zero (level turn) , the vehicle 
depth continues to increase without apparent bounds. This is not a sign 
of instability, but rather an indication that 0 is not a very accurate 
representation of depth rate. 
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VELOCITIES ANGULAR VELOCITIES 

NON-L 1 NEAR SIMULATION Z-PSI NON-LI NEAR SIMULATION Z-PSI 
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NON-L 1 NERR SIMULATION Z - P5 1 NON-LINEAR SIMULATION Z-P5I 







Figure 5-2. Complete nonlinear simulation for \p - z control system 
in vicinity of nominal point. (Cont.) 
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NON-L I NEAR 5 1 MULRT I ON R-ThETR 




TIME 




TIME 




TIME 




Figure 5-3. Nonlinear simulation showing output error and 
response for r - 0 control system in vicinity 
of nominal point. 
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When bearing rates were much greater than the nominal rate, the 
system behaved in a much different fashion. The simulation results in 
Figures 5-4 are for a commanded bearing rate of 2 degrees/s, or about 
double the nominal rate. Note here that the time scale has been expanded 



while the z error shows indications of settling at some time greater than 
400 seconds. Figure 5-5 demonstrates similar results for the r - 0 con- 
troller? the 0 error state has nearly the same shape as the z error state 



complicated and certainly not harmonic in nature. This is no doubt 
caused by the nonlinear behavior of the submarine. 

Ordinarily, performance measures are given in the form of settling 
times, rise times, percent overshoot, etc. These conventional measures 
cannot be applied exactly (unaltered) in this case due to the nonlinear 
behavior that would render such indices meaningless or misleading. For 
instance, a steady-state output error may prolong the measure of system 
settling time, although the output error may have reached this particular 
state in a relatively short time. 

Quantitative comparisons will be made, instead, through perform- 
ance indices that are based on an integral error of the form 



for a step input of y , the commanded output, the integral error becomes 



to 400 seconds. The output error of ip takes almost 120 seconds to settle, 



except for a different sign. Note that the dynamics of z and 0 are rather 




dt 



0 



P.I. 



y 




dt 



(5-1) 
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Figure 5-4. Complete nonlinear simulation for ip - z control system 
away from nominal point. (Cont.) 
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Figure 5-4. Complete nonlinear simulation for ifj - z control system 
away from nominal point. (Cont.) 



NON-L I NEAR SIMULRTION R-THETR 
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Figure 5-5. 



Nonlinear simulation showing output error and response 
for r - 0 control system away from nominal point. 
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This is a particularly useful and convenient measure. As it turns 
out, the change in depth, or depth excursion, in a specified time, t , of 
a submersible that is commanded to maintain a constant depth is equal to 
the depth rate integral error. Equation (5-1) will not be applied pre- 
cisely to the r - 0 controller; however, the resulting depth excursions 
and heading errors provide a common means to compare the performance of 
the two control systems. These results are plotted in Figures 5-6 and 
5-7. Figures 5-6 shows no significant differences between either control 
system's yaw rate integral errors. Several conclusions can be drawn from 
Figure 5-7, however. Both controllers' performance is best near the 
nominal point bearing rate (1 degree/s) , although it can be seen that 
the depth excursion of the r - 0 controller is greater. Again, this is 
not surprising since 0 is only a partial description of the full depth 
expression [ (Eq. (2.14))]. For an equal interval of 1 degree/s on either 
side of the nominal point, the depth excursions of the ip - z controller 
are nearly symmetrical with respect to the nominal point, whereas the 
error of the r - 0 controller is less at lower turning rates. In either 
case, the controllers provide better depth control than otherwise achiev- 
able by manual means, although the rapid increase of depth for commanded 
turning rates greater than 2 degrees/second or less than 0 degree/second 
(opposite turn) may be an indication that the system is becoming less stable. 

It is very unlikely that a submarine would stay in a maneuver like 
this for very long since it would have completed a full (360 degee) turn in 
about 3 minutes. A more realistic operating scenario was simulated. 

Figure 5- 3 shows a series of maneuvers in response to various reference 
inputs listed in Table 5-1. Which are in the form of fast ramps (1-5 
seconds) . This is to simulate the commands an operator may issue through 
the "joy-stick" in maneuvering his vehicle. An interesting observation 
in Figure 5-3 is that commanded changes in 0 have very little effect on r 
as indicated by points 1, 3, and 4; however, the commanded change in r 
does have a substantial effect on the error in 0 , as seen in point 2. 

Point 5 of Figure 5-6 shows the effect of commanding simultaneous changes 
in r and 0 . 
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0.5 0.0 -0.5 -1.0 -1.5 -2.0 

COMMANDED 0 RATE (deg/s) 

Figure 5-6. Integral yaw error for commanded ip 
rate after 100 seconds. 




Figure 5-7. Integral depth error for commanded ^ 
rate after 100 seconds. 
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Figure 5-8. Output error response to various commanded inputs 
for r - 0 control system. 
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Table 5-1. Nonlinear simulation r - 0. Time and the 
corresponding reference inputs. 



Fpirt 



$ 



TIME 

SEC 


r 

(10 deq/s) 


e 

( deg ) 


0.0 


- . 100000E+02 


0. 000000E+00 


60.0 


- • 100000E+02 


0. 000000E+00 


61.0 


- • 100000E+02 


0.200000E+01 


100.0 


- . 100000E+02 


0 . 200000E+01 


101.0 


0 . 000000E+00 


0 . 200000E+01 


200.0 


0 . 000000E >00 


0 . 200000E+01 


201.0 


0 . OOOOOOE+OO 


0. OOOOOOE+OO 


300.0 


0.000000E+00 


0. 000000E+00 


305.0 


0 . 000000E+00 


0.500000E+01 


375.0 


0 . OOOOOOE+OO 


0.500000E+01 


380.0 


0. OOOOOOE+OO 


0.500000E+01 


440.0 


0. 000000E+00 


0.500000E+01 


445.0 


- . 150000E+02 


0 . 200000E+01 


500.0 


- • 150000E+02 


0. 200000E+01 



The full response of the 0 - z compensator to various command in- 
puts are shown in Appendix C. 

5 . 3 Robustness 

To determine how close the controller is to being unstable, recall 

that E (s) was given by 
— pre 1 

E (s) = [G (s) - G (s) ] • G _1 (s) (3-1) 

—pre — — — 

It is clear that for this simulation G(s) above is a function only of the 
system nonlinearities since the hydrodynamic coefficients remain as pre- 
determined constants. The structure of G(s), then, can be found for any 
operating condition. Note in Figure 5-9 that a commanded \ p of -2 degrees/s 
corresponds roughly to 8 degrees of rudder. Since this region of controller 
operation appears to be approaching its performance limit, G(s) will be 
determined by linearization of the submarine dynamics about 8 degrees of 
rudder. 
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Figure 5.9. Steady state rudder required to achieve commanded 
bearing rate. 
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System robustness will then be studied through the inequality 



a [E 
max — pre 



(jw) ] 



a . [I + 
min — 



(G(jw) K ( ju)) ) 




(3-2) 



Figures 5-10 and 5-11 provide a graphical interpretation of Eq. (3-2) above. 

In both instances the maximum singular values of E (s) are equal to, or 

-pre 

slightly larger than the minimum singular values of the inverse return 
difference matrix of the loop-transfer function for the bandwidth up to 
0.1 radians/s. Hence the inequality (3-2) above is not satisfied and 
stability can no longer be guaranteed. 



For the - z compensator of Figure 5-11 a bulge forms in the 

singular values of E (s) . This is most likely due to the presence of 

“pre 

• • 

non-minimum phase zero in the ip - z system since the difference between 

E (s) of either system can only be attributed to the C matrix used for 
—pre ~ 

the outputs. Although the maximum singular value of E^^(s) is not greater 
than the singular values of (I + T ^(s)) near the bulge, it is typically 
at these higher frequencies and above that other modeling errors (such as 

unmodeled dynamics) are likely to occur in an actual submarine. The combina- 

• # 

tion of the unmodeled errors in G(s) and the bulge in E (s) for the ip - z 

— —pre 

controller is more likely to result in an unstable system than, say, rela- 
tively flat singular values of E^ re that appears in Figure 5-10. 



5 . 4 Gain Scheduling 

Gain scheduling involves the use of more than one set of matched 
Kalman filter gain matrices (H) and control gain matrices (G) to provide 
the compensator with the ability to extend acceptable performance and ro- 
bustness requirements over the widest possible range of operating condi- 
tions. These operating conditions for the rate controller design can be 
expressed here in 3-dimensional "volumes". This concept is shown if 
Figure 5-12. The specification of speed, rudder deflection, and stern 
plane defection define a single point. On this diagram, one can plot 
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SINGULAR VALUES (LOG MAGNITUDE) 



Ep re , (I + T ' 1 ) FOR R-THETR SYSTEM 




Figure 5-10. Inverse return difference of recovered transfer function and 
premultiplicative error; r - 9 system; singular values vs 
frequency. 
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SINGULAR VALUES (LOG MAGNITUDE) 



I 



c 



L pre 



+ T' 1 ) POR PS I - Z OCX COMPENSATOR 




* i 1 — r -rr ? t — r ' r~rn * — T — — i — t — i r n i ; 1 1 — t — r i m 1 i i — ! — i — r~ rrn 

10' 3 iO' 3 10' 1 10° io 1 io 3 

FREQUENCY (RAD/SEC) 



Figure 5-11. Inverse return difference of recovered transfer function 

• • 

and premultiplicative error; ip - z system; singular 
values vs frequency. 
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Figure 5.12. Schematic showing operating volume for determination of model 
robustness and acceptable performance near a nominal point. 
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regions where, about a nominal point (s), the multiplicative error can be 
determined and robustness checked in a test similar to Eq. (3-2) above. 

This is called the region of guaranteed stability. Another volume on the 
same plot is reserved for performance. It includes all points about the 
nominal point that result in acceptable performance. The intersection of 
the stability region and the performance region defines the outer bounds of 
the operating range of the particular set of gains. Many nominal points 
together would be required to cover the full possible operating range of 
the submersible. 

Gain scheduling can then be accomplished by two possible approaches: 

(1) Establish large overlapping volumes around nominal points 
where particular set of gains (including corresponding models 
of the A, B, and £ matrices) are satisfactory with respect 

to performance and robustness. The shift of gains from one 
region to another would be through a band-bang implementation. 
This procedure, of course, requires implementation schemes 
that would allow the gain shift to occur without violent 
transients of the control surfaces. 

(2) Select nominal points to span the volume, and schedule the 
gains of the compensator through a quadratic or least squares 
fit of the operating conditions through the nominal points. 
Lively [14] has successfully implemented such a methodology. 

It becomes difficult, however, to evaluate the robustness by 
methods used in this thesis. One must somehow ensure that 
all quadratically obtained gains result in enough robustness 
to overcome the unmodeled dynamics. 

5 . 5 Other Notes 

When q, the control gain weighting index is increased over the value 
of 10 used in this design, the gain entries of £, of course, increase. The 
result is that: 



99 



(1) The bandwidth of T(s) increases slowly to match the crossover 
° f — KF* 

(2) The control surface deflections for maneuvers similar to those 
performed in this chapter increased, and tended toward satura- 
tion as q increased further. 

(3) The maximum rate limits of the actuators were exceeded. 

(4) There was no noticeable increase in performance. Therefore, 
the value of q was kept as small as possible consistent with 
meeting the minimum bandwidth requirements. 

5 . 6 Chapter Summary 

The performance comparisons of and between the two compensator 
designs are far from complete. It was the intent of this chapter to show 
how one could begin to evaluate a design; a full evaluation of each con- 
troller would require many more calculations to be completed. One might 
conclude at this point that linearization required to obtain the C matrix 

used in the ip - z controller has resulted in a somewhat less robust design. 

• • 

However, the \p - z controller provides more precise control of turning and 
depth rates which is the ultimate control objective of this design. 

The question of robustness away from the nominal point chosen in 
this design was only partially ansered. Nonlinearities most likely pro- 
vide the greatest source of modeling errors. Other dynamics that arise 
from errors in the hydrodynamic coefficients and unmodeled dynamics are 
not known at this time but are required in order to complete the robustness 
picture before a practical design can be produced for a real submersible. 

The control gain matrices for both control systems are found in 
Appendix B. 
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CHAPTER 6 



CONCLUSIONS 



6. 1 Summary 

Multiple -input multiple output control system design methodologies 
has been successfully employed to provide a means for an operator of a sub- 
mersible to control heading and depth for complicated maneuvers. This 
study is far from a completed design. However, results thus far are en- 
couraging. Of significance this thesis has: 

(1) Provided one of few available examples of a complex LQG/LTR 
Model Based Compensator designs evaluated on a nonlinear 
simulation. 



(2) Developed a loop- shaping technique, using integrator aug- 
mentation in the command input channel, in which singular 
values can be brought arbitrarily close together and then 
recovered through natural scaling of the output variables. 

(3) Compared the performance and robustness of rate control 

designs based on a constant C matrix (the r - 0 controller) 

• • • • 
and a linearized C matrix (the \p - z controller) . The ip - z 

controller, which directly controlled the heading rate, \p , 
and the depth rate, z, provided better control of the sub- 
marine depth in a turn than the r - 0 controller design. 

Heading rate control was essentially equal in both designs. 

• • 

However, the linearized C matrix of the \p - z design showed 
potential for less stability as the vehicle states deviated 
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from the nominal point. The non-minimum phase zero of the 

• * 

]p - z controller was outside the desired bandwidth and did 
not effect the loop-transfer recovery process. 



6.2 Recommendations 

Much more work can be done to further provide practical uses of 
the MIMO LQG-LTR methods used in this thesis and to take a deeper look 
into achievable submarine performance through rate control systems. Some 
topics not thoroughly covered in this thesis that could provide a basis 
for further research include: 

(1) A demonstration of stability and performance as a function of 
state variable scaling to see if the natural scaling pro- 
cedures developed in this thesis effect the actual performance 
bandwidth over unsealed state variables or designs with differ- 
ent scaling methods. 

(2) Determining an error norm that combines both the known errors 
due to the model deviations from linearity as well as errors 
that estimate the tolerances of the hydrodynamic coefficients 
and other modeling errors in order to provide a better picture 
of the achievable compensator robustness. 

(3) Checking the performance of the model in response to environ- 
mental noise. 

(4) Developing a practical method of gain scheduling so the con- 
troller provides good performance in all attitudes and speeds. 

(5) Determining how scheduling the sail planes (the k factor in 
Eq. (2-11)) effects the performance and stability of the sub- 
marine . 

(6) Applying the design procedure used in this thesis to develop 
control systems for different fin configurations, to possibly 
include "Y" and "X" sterns, and differential control of a 
cruciform stern. 
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Table A-l. Dynamical response terms (Abkowitz) . 

Body Properties 

Length 

Geometry 

Mass 

Moments of Inertia 



Properties of Motion 
Orientation 

x , y , z fixed coordinate system, reference at surface 

o o o 

<J> , 0, ip angular rotations about the x, y, and z axis, 

respectively 

Body Motion 
~y 

U linear velocity vector 

ft angular velocity 

U linear acceleration vector 

♦ 

ft angular acceleration 

5 control surface deflection 

6 control surface velocity 

6 control surface acceleration 

n propeller angular velocity 

n propeller angular acceleration 



Fluid Properties 

p mass density of fluid 

y viscosity 

t surface tension 



g acceleration of fluid 

p pressure of fluid 

Py vapor pressure of fluid 

E elasticity of fluid 
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A. 2 



Nonlinear Equations of Motion 



Explanation of Conventions used in Hydrodynamic Coefficients 

Hydrodynamic coefficients used in the nonlinear and linear equations 
of motion are generally written in the form 

X , 
ab ... 

where X represents a force or moment and the subscripts a, b, ... represent 
an angular velocity, linear velocity, rudder deflection, etc. The con- 
ventional meaning of coefficients written in this form is 



which are taken to be Taylor expansion coefficients. 

Moments of inertia are written in similar form, but all subscripts 
following retain the usual meaning of coordinate system. 



X 



ab. . . 
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q.|<n q.|<n 



1 



SURGE 



m[u + qw - rv - x Q (q 2 + r 2 ) + y Q (pq - r ) + z Q (pr + q) ] * 

p 4 r 2 2 i 

^ Z [X q + X r + X rp 
2 L qq rr rp J 

+ Z 3 [x*u + X vr + X wql + Z 2 [x u 2 + X v 2 + X w 
2 L u vr wq J 2 L uu w ww 



i 2 U 2 [x.„ 6r 2 + x. . 6s 2 + X....6b 2 ] + (PROPTH - DRAG) 



2 2 2 2 2 2 2 
z (n-i) [x v + X w + X u ilk + X. u 6s 1 
L wo wwrj 6r5rri *<**<*r» -J 



6s6sr) 



- WTOT sin 0 



2 2 

PROPTH - DRAG = u f c . n + b . n + a . ) 

^ 1 i 
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2 



SWAY 



m[v - pw + ur - y G (r 2 + p 2 ) + Z( ,(qr - p) + x Q (qp + r)] = 

j 4 4 [Y«r + Y*p + Y i ip|p| + Y m pq + Y^qr + Y_ sign (r,v) /v 2 + w 2 ] 



pq qr" 



+ I ^ V*P + Y r Ur + V P + Y |r|6r U | r | 6 ^ + 

7 ^[y^u 2 + Y „ uv + Y | r |v v/v 2 + w 2 +• Y. u 2 6r + 



4* Y vw 
vw 



,/2 2+ 2 i 

u// u 4* v W j 4* WTOT sin <J> cos 9 4 * 



+ (n-1) i 3 Y ur + Z 2 Y uv + £• £ 2 Y. u 2 6r + 

7 ™ 7 vt) 2 5 rn 



r n 



^ Z^Y | | Wv4-w^ } 4- Y 4- Y 

2 v r Tl J CRFL VTX 



Y CRFL : ^ 2 (~ l 



P r 2 ZWW 



f D(x)dx 



)} /» H(x)v(x) /v 2 (x) 4* w^(x) dx 



j>_ „ 2 Z w 

"VTX ’ 1 2 



* r P 2 

Y : { — Z 

VTY l 9 



i r f w — 

} J w(x)v fw (t - t ) dx ; not included in 

IBM program. 



(x_ - X ) X 

fw so vs 



sign(r,v) = r • y^j- 



t y : time it took vehicle to travel a distance x in the axial direction 
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3 



HEAVE 



• 2 2 • 

n[w + vp - uq - z G (p + q ) + x Q (rp - q) + y Q (rp + p)] 



rr £ 4 [ Z*q + (z r 2 + Z rp) 
2 L q v rr rp CJ 



/2 2 2 
/u + v + w 



q- Jl 3 [z.w + Z % 
2 L w vr 



vr 



u „ w / 2 2 | | 

7==== * «| q l FT /v lql * 

/U + V + W 



+ Z g uq + Z| g | u | q | 6 s] + 



P „ 2 r U „ 2 _ VU 

+ ^ Z z 

2 L vp 



/ 2 2 2 
/u + V + w 



vp + Z^u + Z 



w 



/ 2 2 2 
/u + V + w 



/ 2 2 2 2 || 

+ Z t | /v + w W + Z uw + Z. u 5s + Z p ,u 5b + Zi iu w + 

w 5s 5b I w | 



+ Z^lwl /v 2 + w 2 ] + WTOT cos $ cos 9 + 



r p 3 p 2 p 2_ / 2 2 

+ (Tl-1){^£ z uq + ^ Z Z uw + — £ Z i i w/v +w + 

1 2 qn 2 wn 2 w w D 



+ 2 1 Z 5sn U 5s ^ + Z CRFL * Z VTX 



CRFL 



= - 2- f C D(x )w( x ) /v 2 ( x ) + w 2 (x) dx 

s * H 



2 i d 



VTX 



fw 

2- £ c f v(x) v- (t - t )dx ; not included in the IBM 
2 y J fw x 

J x program. 
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4. ROLLING TORQUE (<j>) 



X x P + (I z ‘ X y )qr ‘ { ' r + w)I xz + (r2 " q2)I yz + (pr ' q)I xy + 



m[y_(w - uq + vp) - z (v - wp + ur)] 
G G 



T ^ [ K *P + V + K nlo|Pkl + + K rwN] + 



P P 



qr pq 



+ % Z 4 f K* v + K ur + K upl + 

2 L v r P J 



+ 5. [ k . u 2 + K uv + K | j v/v 2 + w 2 + K_ u 2 6r + K vw + 

2 L * v v j v j 5r vw 



+ (rj -1 ) K u ] 

*T1 

- (z WTOT - ZB) cos 9 sin + (y WTOT - y B) cos 9 cos <j> + 
G B G B 



+ [u^ + V 2 + w 2 ] [k^B^ Sin 4<j> t + sin ftfrj + 



2 



fw 



+ — ■ l C Z / w(x) V (t - T )dx 

2 v J w x 



vs 



v = v + x r 
t t 



w = w - x q 
t 



= tan' 1 (- 



-1 72 2 

tan /v + w^/u 
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5 



PITCHING TORQUE (9) 



I q + (I - I )rp - (p + qr ) I + (p 2 - r 2 )I + (qp - r)I + 
y x z ^ xy zx yz 



[z G (u - vr + wq) - x Q (w - uq + vp)] = 
2 u -i 

J + 



= Z^ f M*q + M r 
2 1 q rr 



/ 2 2 2 
/u + V + w 



+ | Z 4 [M.w + M| w | q q/v 2 + + M q uq + N | q | & g U | q | 6s ] + 



2 2 



+ K [M,u 2 + M 5s v 2 



f 2 2 

+ M | jw/v + w + M | iu w 

/ 2 2 2 w w w 

/u + v + w 



^ ^ / 2 2 

+ M^u^Ss + M^^u^5b +M | W | u | w | +M W | W | | w | + w ] 



- (z WTOT - z B) sin 9 - (x WTOT - x B) cos 9 sin 
G B G B 



+ (n-1 ) [ £ 4 M uq + Z 2 M uw + ^ i 3 M i i w/v 2 +w 2 + M. u 2 5sl 

1 2 qn 2 wr\ 2 w w n 5sn J 



+ M + M 

CRFL VTX 



M CRFL * 2* / C d X ° (X) W(X) + w2 < x > dx 

Z 



VTX 



A- 

fw 

- rr Z C / x v(x) vf (t - t ) dx ; not included in IBM, 
2 v J w x 

x 

vs 



1 12 



6 



YAWING TORQUE (iji) 



V + (I y " r x )pq ~ (q + rp)I yz + (q2 " p2)I xy + (rp ” P)I 2X + 

+ m[x_(v - wp + ur) - y_ (u - vr + wq)l = 
u G G J 

'x£^fN»r + N«p + N qr + Nr i i r I r I + N pq] + 

2 L r p qr r II pq J 



+ 2-£^rNii r/v 2 + w 2 + N wp + N ur + N up + N*v + Ni i. ulrl 6r] + 

2 L | v| r wp r p v |r|5r I I J 

] - 



[N u 2 + N i | v/v 2 + w 2 + N 
9 L * v V I 



UV + N 
V VW 



/ 2 2 : 
/ U + V + w 



+ (n -1 ) £^N ur + £ 2 N uv + ^ ^ N i i v/v 2 + w 2 + 

l 2 m 2 vti 2 v v n 



* s *\ a ^ 



+ (x g WTOT - XgB) cos 9 sin <J> + (y G WTOT - y 0 B) sin 9 + N CRFL + W VTX 



CRFL 



|c / X H(X) V ( X ) /v 2 (x) + w‘ 



(x) dx 



VTX 



fw 

- l C / x w(x) vf (t - t ) dx ; not included in IBM 
2 v ; w x 

X 

vs 
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N • 



KINEMATIC RELATIONS 



2 2 2 2 
U = u + v + w 



-u sin 9 + v cos 9 sin d> + w cos 9 



p + iji sin 9 



q - ip cos 9 sin 
cos <p 



♦ = 



r + 9 sin <p 
cos 9 cos <p 



cos <p 
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A. 3 



Parametric Linearization 



I. AXIAL FORCE (SURGE) 



[m - Xu] Au - my.Ar + mz Aq = 
G G 



Au [ 2u fx + a. + b,n + c.n + X_ 5r + x 
L o L uu i i i 5r6r 5s6s 



2 V X 2 

X- p 6s + 



X 5b6b 5b2 + (n o- 1)(x 6r6rn 5r2 + X 6s6sn 6s2)} l + 



+ (n -Dx } v + 

O WTV O 


(X + 


m )r o ] + 


vr 


+ (n -Dx }w + 
o wwn J o 


(X 


nOqJ + 


wq 


- (mz -X )r 1 
) G rp o J 


+ 




» m)w - my d + 2(mx„ + 
O Go G 


X )q 1 + 

qq o J 



+ Ar [ (X + m)v + (X - mz )p + 2(mx + X )r ] 
L vr o rp G o G rr o J 

- A0 [ (WTOT - B) cos 0 ] 



+ A5s [ 2( X + X (n -1))6s u 2 ] 
L 6s6s 6s6sn o o o J 



+ A5b [ 2X u 5b 1 
L 5b6b o o J 



2 

+ A5r [ 2( X x x + X x . (n -1 ) ) 5r u 1 
L 6r6r 6r5rn o o o J 



An [u (b. + 2c.n ) + X v + X w 
1 o i i o wn < 



2 + x w 2 
wn o wwn o 



+ (x 5s 2 + X 5r 2 )u 2 ] 
^ 5s5sn o 5r6r o' o J 
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II. LATERAL FORCE 



[m - Y* ]Av - (mz_ + Y» ) Ap + (mx_ - Y» )Ar 
v G p G r 



Au f 2u { Y + V. fir + (n -1 )Y. fir - ^ 

*• o *• * fi-r n rtrn 7 



y u v w 
vw o o o 



6r 'o 5rn 2 r 2 2 2> 3/2 

u + v + w ] ' 

K o o o' 



v f Y + (n -1 )Y + 

O L V O VT| 



y w 
vw o 



/ 2 2 2 

/ U + V + w 
o o o 



} + + 



P O 



+ r 



{ (Y r - m) + Y 



r fir r 



fir + (n -1 )Y 1 + 

o rn J 



y w 

r r vw o r 

+ Av |u {Y + fl - 

L 0 1 V , v 

/ 2 2 2 

/u + v + w 
o o o 



2 2 2 

u + V + w 
o o o 



) + (n -i )y + 

J o vn 



V {. 

O 



[2 T 

/ v + W 
o o 



(y [ |+ (n-i) Y | i )•( i + + - )v 

v v|v| V V n ; v v 2 7™2 2 ° 

o /v + w 

o o 



a Y 



+ Y q + 



CRFW 



vq o 3 v 



+ Aw [u { 



Y v 
vw o 



/ 2 2 2 

/u + V + w 
o o o 



(1 



w 



2 2 2 

u + V + w 
o o o 



)} 



w r 

v { ( Y I I [■■■— ■ + Y I I + (n -1)Y I I 

° ri — 2 v l r l K\ v M ° v M n 

/V + w 
o o 



3 Y 



+ p { Y + ml + Y r +pm + 
o l wp J wr o o 



CRFW 



3 w 
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+ Ap fu Y +w(Y +m)+p ( 2my _ + 2 — — Y j i) + q (Y - mx ) 1 + 
c L o p o wp G P Q p|pr o pq 



pq 



3 Y 



Aq f Y v + p ( Y - mx ) + r (Y - mz ) + 

^ L n nn ft r> rrr ft 



CRFW 



i L 

+ Ar [u { (Y - m) + 2Yi i t ° | 5r + Y (n -1)} + v {y i i I — 
L o l r r 5r To rn o 1 o l v r |v 



l/2 2 

• — / v + w } + Y w + q (Y - mz _ ) + r • 2m v + 

r O O J wr o nr ft o ^ ft 

O 



3 Y 



CRFW 



wr o o qr 

+ A 4> [ (WTOT - B) cos 9^ cos <j) Q ] 

+ A9 f - (WTOT - B) sin 9 sin <j> 1 

r» n J 



3 r 



2 

+ A5r [Yi I u r + {y, + (ri — 1 ) Y _ }u ] 

L r 5r 0| ol 1 5r o 6rn J o J 



2 2 2 

+ An [ Y ur + Y uv + Y i iv/v + w + Y p u 5r 1 
L rn o o vn o o r v n o o o 5rn o o J 



where: 



3 Y 



CRFW 



3 v 



p r „ 2 ww > r , , ' < 

= ” f-A r —- / — rr J /* h x — 

2 ^ J D(x) dx ; J i 



(w 2 (x) + 2 v 2 (x)) 
v o o ' 



dx 



/v 2 (x) + w 2 (x) 
o o 



3 Y 



CRFW 



3 w 



- * v (x)w (x) 

0,2 ww > r . x o o 

f .< i 7 _- — . — — ) f H( x ) dx 

2 v J.D(x) dx ; J £ r 

/v (x) + w (x) 



✓v (x> 
o 



3 Y 



CRFW 



3 q 



2. (i 2 , 222 ) f x 

2 L / D(x) dx j J 2. 



H( X ) V ( X ) w (x) 

o o 



/v 2 (x ) + W 2 (x) 
o o 



dx 



3 Y 



CRFW 



3 r 



!(-‘ 2 T 75 ^> * ■'*' 



( 2v 2 (x) + w 2 (x)) 
v o o J 



J7 2 



(x) + w (x) 
o o 
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III. Normal Force (Heave) 



(m - 



Z* )Aw + my_Ap - (Z« + mx )Aq 
w G q G 



Au[2u 0 { z 5s 5s + Z 6b 6b + Z, + (n n - DZ^ n 6s} + 



5sn 



+ v { (Z v + Z _p_ + Z r ) f 
0 1 w 0 vp2 0 vr o ^ 



/ 2 2 2 (u? + v 2 + w 2 ) 3/2 

/u Q + v 0 + w 0 V o 0 0 J 



•)} 



I 



+ w f z + z | I + (n - i )z } + 

0 l w |w| 0 wiy 



+ pJ (Z p^ + Z r.) f 
o 1 ppO rp 0 v 



/ — = r e 2 2 20/2 

/ 2 2 2 u + v + w 1 7 

/ u + v + ^ 0 0 0' 

0 0 0 



)} 



+ q.fz + m + Zi 

0 l q q q. 



63 + (r >n " 1)Z .n} 



qn J 



+ r_{z r _ f 
0 l rr 0 1 



/— ( 2 2 20/2 

/ 2 2 2 (u + v + w 1 

/ u + v + v/ v 0 0 O' 

0 0 0 



■)}] 



Av [ u 0 ( 



(Z r + Z ^p + 2Z v ) 
vr 0 vp2Q w 0 

f ~ 2 2 2 

/ U + V + 

0 0 0 



2 2 2 
v f Z r + Z p^ + Z r p^ + Z v r_ + Z v + Z v 1 
0^ rr 0 pp 0 vp 0 0 vr 0 0 vp2 0 0 w O' 

( 2 2 20/2 

( u + v + w 1 
^ 0 0 0 J 
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III. Normal Force (Cont. ) 



v ° {Zw M KT 



,2 2 / 2 2 
/ V _ + 

0 0 0 0 



(Z 



I w + Z I |W A + (n - 1)Z I I w ) } + p fz 
I 0 1 w w 0 0 w w n 0 J 0^ vp 3r 

w Q 

+ Aw[u_{z + Z | i -| 1 - + (ti - 1 ) Z 

L 0 l w w W 0 wri 



CRFW 



)] 






( 2 2 2 ^ 3/2 

[u^ + v + w J 
^ 0 0 O' 

2 2 2 
fz + Z r + Z rp rt + Z vr + Z _ v_p_ + z r 1} + 
^ pp 0 rr 0 rp 00 vr 0 0 vp2 00 w 0' J 



2 2 

+ v ~{ z I I /I + + Z 

0 1 w w * v n' n ' ^ 



w v 
0 0 



0 0' ww 



i i J 2 

| w o| /v o + w 

t i u i i u | i u i 0 

+ w 0 l z w| q | — * 

| q | , 2 2 ,2 2 

/ V + i + ** 



+ (T1 o ' n 



w w m o 



j 2 2 

0 0 



q_ . 

1 0| | H 0 I 0 | 



0 0 



0 0 



3 Z 



r __ ww _ , > CRFW -i 

( 2Z | t w + -i j- * 2w (n - 1 )Z | | • 2w J + + 

^ w w 0 w 0 o w w m 0' 3w J 



Ap[u 0 



,2 2 2 

/u + v + 

0 0 0 



( 2Z p rt + Z r_ + Z v ) + 

pp 0 rp 0 vp2 0 



+ V 0 (Z vp" m) + P 0 2mZ G “ q 0 my G ' r 0 mX Gl + 
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III. Normal Force (Cont. ) 



q o 

+ Aqf u { Z + m + Z i i. -j 1 - 6 s + (n ^ - 1 )Z } 

L 0 l q q 6s q 'o qn J 



pi 2 

/ v + w _ 

r _ 0 0 0 1 

+ w n Z j j j 1 1 f- | + q 2mz - r_my + 

0 w | q | ro| | q o| 0 G ° g 



3Z 



Ar [ u n( 



1 



0 L 7 



, 2 2 2 

/u + + V/ 

0 0 0 



( 2Z r + Z p + Z v ) } 
rr 0 rp 0 vr 0 J 



3Z 



CRFW 



- p^itdc - q^my^ + 

*0 G ^0 *G 9r 



] 



+ A cf> [ -( WTOT - B) cos 8 q sin 4> Q ] 

+ A9 [-(WTOT - B) sin 0 Q cos <J> Q ] 

+ A6S t (Z 6s + ( % - 1>Z 6sn )U 0 + Z |q|6s U o| q o|3 + 

+ A6b t Z 6 b u o} + 



/ 2 2 2 

+ An I" Z u (J + Z u w^ + Z I I w Vv + + Z. u^6s 1 

L qn o o wn o o w w n o o o 6sn o o J 



where : 



3Z 



CRFW 



a v 



- o’ (”* 2 



ww 



2 ^ j l D(x) dx- 1 'Jt 



w (x)v (x) 

•) /, fix) —2 2 a* 



/v^(x) + W Q (X > 



CRFW 

3q 



3Z 



CRFW 



a w 



p f n 2 WW >1 r rv/ . 0 

2 J D(x) dx^ ^ i ° (X) - 



(v Z (x) + 2w Z (x)) 



dx 



/vq(x) + w Z (x) 
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IV. 



3z crfw 

3q 



* 





a$ U x D(x) 



(v 3 (x) + 2 Wq(x) ) 

■ dx 

/v 2 (x) + w 2 (x) 



3 Z 



* 



CRFW 



3 r 




ww 



/^D(x) dx 



) /, 



D(x) w Q (x) v Q (x) 

x — dx 

/Vg(X) + W Z (X) 



Rolling Torque 



-[mz^ + K» ]Av + my_Aw + [I K« ]Ap - I Aq - [I + K« ]Ar 
G v G xp xy xzr 



Au[ 2u { K* + K. 6r + (n n - 1)K. + 

L 0 l 6r 0 0 *n 



+ S tO (K 4t 3in 4 ^n + K a*- sin + 



to 8t 



to J 



2 w 2 

+ V o{ K v - 28 to J-T + “T * (K 4t Sin 4 *t 0 + K st Sin 
v o v o 



+ p K + q my„ + r^mz^ + M Z_ + K ) r 
*0 p ^0 J G 0 G o G r o 



Av[u {- 



2v 8 

o p to 



/ 2 2 

v to + w to 



u o 8 to w to 

(K 4t Sin 4 *t0 + K 8t Sln 8 ^t0 ) + “2 T 

v to + w to 



( 4K sin 4<j> + 8K 0 ^ sin <j> J + K } + vJk i if 

4t ^t0 8t Y t0 v 0 L v v 



/ ^ 

Ji + _ + + 



v , 2 2 

0 /v Q + w Q 



+ 8 tO^ (K 4t Sin 4 *t0 + V Sin 8 ^tO )2v tO + 
+ W t0 (4K 4t COS ^tO + K 8t 8 COS + 

+ Vuv - P 0 my G + K vq q 0 l + 
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IV. Rolling Torque (Cont. ) 



^to^to , 

+ Aw|uJ (K sin 4<j> + K sin 8<f> ) 



Aw[u Q {- 



2 2 

v , + *7 

to to 



u o 8 to v to 



2 2 v 4t 

v _ + w „ 
to to 



(4K 4fc cos 4<j> tQ + 8K 8t cos 8<0^ n )} 



to' 



+ V,JK 

o l ' 



K I I 
v v 0 



vw 



/ 2 L 2 
/v A + v/ 

0 0 



+ 8 to{ 2W tO (K 4t * Sln 4 **n + Sin 8< ^n> 



to 8t 



to 



“ V tO UK 4t C0S 4 *t0 + 8K 8t COS + 

+ V K wp * "V + VwJ + 



+ Ap[ u k - v my + w{k - m z\ + 2k i i Ip . 

L 0 p 0 G 0 l wp G J p p | 0 1 + 



+ q ( I + K ) - r I ] 
0 yz pq 0 xy J 



r f , ^toWtO ^ „ 2 , 

+ Acj [ u nl my n - ( ~ =r + 



/ 2 2 

✓v „ + w . 
to to 



(K 4t sin 4^ 0 + K 8t sin 8^ Q ) + 

4u v x 

+ ------ 1 " (lc 4t COS 4 *t0 + K 8t 2 COS ^to^ + 

v t0 + W t0 

* Wl * •’ol 1 ,. * W - W * r oI a y - V * 

* K qr * 46 tO v tO x tO <,C 4t “ S 4 *to * V 2 c °* 8 *to’! * 
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IV. Rolling Torque (Cont. ) 



Ar [ u 0 { mZ G + K r + 



26 to v to x to 



/ 2 2 

/v „ + w 
to to 



(k 4t Sin 4 *t0 + k 8t Sin ^tO 5 + 



U 0 S t0 W to X to 

+ ' " (4K 4t *t0 + 8K 8t COS 8 *t0^ + W 0 K wr ' Vxy + 

v to + w to 



q_{ (I - I ) + K } - r 21 + 

0 L y z qr J 0 yz 



+ 2 v t0 X t0 6 t0 (K 8t Sin ** to + K 8 t Sin ^to 5 

+ 6 tO W tO X tO (4K 4t COS 4 ^t0 + 8K 8t C ° S ^to’l 
+ A4> [-(y Q WTOT - y g B) cos 0 Q sin 4> 0 - (z^WTOT - Zg B ) cos cos 9^] + 

+ A 9 [ (z^WTOT - z^B) sin 9 Q sin 4> Q - (y G WTOT - y B B) sin 0 Q cos 4> Q ] + 

* Wr I V*o) 

* M [%»Sl 
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V* Pitching Torque 



mz^Au 



(mx^ + M* )Aw 
G w 



I Ap + [I - M» ] Aq - I Ar 
xy y q yz 



= Au[2u {M* + M 5 s 5s 0 + m 5b 6b 0 + (n 0 - 1 >M 6s} + 



2 2 

+ v* 

0 0 



+ V 0^ M vr r 0 + M vp P 0^ + M w V 0^ [2 2 2 ] 3^2 

K + v o + V 



+ w {M + Mi I + (n - DM } + 

0 l W I w I W Q 0 wn J 



2 2 

V + V/ 

0 0 



+ P 0 (M pp P 0 + M r P 2 r 0 ) [“5 2 2^2 

( u 0 + v 0 + V 



q J M - Hoc + M 

^ n*- rr C, 



6s 



1 0 1 * *q G ' |q 6s q 0 0 qn J 



n + (n A - DM } + 

0 0 qn 1 



M r 2 

r M r v + v/ 

0 rr 0^ 0 0 J 

r 2 2 2>3/2' 

1 + V + w 1 

^ 0 0 O' 



h 



Av [ 



{m ( 2v (u^ + w^) + v^) 
1 w v 0^ 0 0 J 0 J 



f 2 2 2>k 3/2 

[u. + + w J ' 

^ 0 0 0 J 

2 2 

- v f M r + M p^ + M r^v^ + M v,d + M d 1 } 

0^ rr 0 pp^O vr 0 0 vp CrO rp2 Cr0 ;j 



2 2 2 

/u^ + + v/ 

0 0 0 



(M r _ + M p + M r ) + 
rr 0 vp O vr 0 



12 2 

V v + v/ 

0 0 



( M 1 1 q_ -f- M I j w_ + M |w I + 

1 w q 0 w w 0 ww| Oj 



(n « - DM w 1 + p_i 
0 wwr) 0 J 0 



3M, 



ntx + r mz,_ + 



CRFW 



_ T i. _ ilU- _ T a 

GOG 3v 
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V. Pitching Torque (Cont. ) 



+ Aw[u 

( 

{ 

H 

+ Ap [ 
+ Aq[' 

+ Ar 



- w„ 



| U 

L 0 V 2 2 2a 3/2 

^ 0 0 0 J 

2 2 2 

M p + M r + M vr + M v d + M v,d + M v 1 
pp 0 rr 0 vr 0 0 vp 00 rp2 00 w 0 J 



■ M + Mi j f t + (n - 1 ) M } + 
w w w I 0 wn J 



[2 2 
/v o + w o 



Mi | q„ + M 



r ^ 2 2> 

( 2w 0 + v Q ) 



w I q 0 w w w 



+ M 



r 2 2 a 

O 0 + V 0^ 



ww I w . 



3 M^ 



h (n - Dm 

‘o wiwin 



f ^ 2 2 a 

(*0 + V, ~“CRFW 

T 0 1 - 1112 G q o + — a^r 



/ 2 2 2 

/ + w 

0 0 0 



( 2M p + M v + M *r_) + mx v 
pp 0 vp 0 rp2 0 GO 



2p I - q I + (I - I + M )r ] + 
0 zx 0 zy z x rp 0 J 



r q 0 r i 

i_ f M - mx_ + (ri - 1 )M + Mi i x *+-- — t- 6s } 
0 l q G 0 qn q 6s q_ 0 J 



Z~*\ 



f v_M | I /l + — - mz w - I p. + 2M II q + 
0 j w j q ^ 2 GO yz 0 q|q| I 0| 



9 M 



fir 



CRFW 



xy 0 3 q 



,2 2 2 

/u o + v o + w o 



(M v + M n + 2M r ) + mz v + 
vr 0 rp2 0 rr 0 GO 



3M 



CRFW 



+■ p (I - I + M ) + I q. + 21 r. + r— 

*0 z x rp xy 0 xz 0 3r 
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V. Pitching Torque (Cont. ) 



Where 



+ A<J> [-(x^WTOT - x b b ) cos 9 q sin $ Q ] + 

+ A9 [ (x W - x B) sin 9 cos <f> - (z WTOT - z B) cos 9 1 + 
Cjb UUG B 0 J 

+ A6s t (M 5s + M 6sn (n 0 " 1))u 0 + M |q|6s U o| q oN + 

+ A6b C M (Sb U 0 ^ 



* 3 m 



CRFW 

a v 



+ An t M qn u o q o + M wn u o w o + M w|w|n w o /v o + w o + 

+ M 6 sn u o 6s ol 

= e. r z 2 .—jsl ) r 

2 ^ J D(x) dx J J l 



D( x ) w Q (x) v Q (x) 
x dx 



/Vq(x) + w q(x) 



a m 



CRFW 

3w 



K' 2 J 7 ( 7 h 7 ) 



(vq(x) + 2Wg(x)) 
/vq(x) + Wg(x) 



dx 



a m 



CRFW 

3q 



p ( n 2 ww > f 2 x 

2 ^ J^D(x) dx^ U X ° X 



( v 0 (x) + 

/vq(x) + w^(x) 



3 M 



CRFW 
3 r 



p 2 ww > r 2^ /x 0 

2 v J D(x) dx ; J Jl 7-r 



w^(x) v Q (x) 



dx 



/v^(x) + w q(x) 



126 



VI. Yawing Torque 



-my^Au + (mx_ - N* )Av - (I + N* )Ap - I Aq + (I - N* )Ar 
G Gv zx p yz ^ zr 



Au[2u 0 {N* + (N 6r + (n 0 - 1)N 5 )6r 0 } + 



rn 



( 2 2^1 
w [ v + w 1 
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^ 0 0 Q J 



+ N p + r_[*N - nut + Ni | _ — 

p^O 0 L r G r 6r r 



6r + (n A - 1 )N 1 

o o m J 



( 2 2^1 

w fw + u 1 
0 O' 



+ Av u f N + N - r — rr ~ + (n - 1 )N } 

L 0 l v vw r 2 2 2^3/2 0 vtv 

fu^ + v + w 1 ' 
v 0 0 0' 

{ N | | r+2N j | v_ + 2(ti _ - 1 )N i | v 1 + 

1 v r 0 v v 0 0 v v n 0 j 



/ 2 2 

/ v + w^ 
0 0 



/ 2 2 

/ v + w^ 
0 0 



f N I iw rt + (n - 1 )N | I w} + 
1 v v 0 0 v v n 0 J 



3N 



+ N q n - my r + 



CRFW 



vqO GO 3v 



Awl* u N 
L 0 \ 



(2 

v o^ u o + v o^ 

vw , 2 2 2> 3/2 

I + V,+ W 1 ' 

v 0 0 0 ; 



/ 2 2 

/v rt + W^ 
0 0 



In I | w + (n - 1 )N I I w } + 
1 v v 0 0 v| v\T) O j 

_(n 



N | | w 

v r 0 



+ r J N + 

wr 



3 N 



+ mx d_ + N p_ + my q + 



CRFW 



T 2 2 

v v + w 
0 0 



G^O wp^O 1 G^O 3 w 
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VI. Yawing Torque (Cont. ) 



where: 



+ Apf N u + w ( mx + N ) + 2p^I + q • 
p 0 0 G wp xy 4 0 

(N + I - I + mv + N ) + I r 1 + 
pq x y J G pq yz 0 J 

+ Aq[ N v n + my w n + (N + I - I )p^ 

L vq 0 J G 0 pq x y *0 



-21 q rt + r (N - I ) + 
xy 0 qr zx ; 

r „ 



3 N 



CRFW 



aq 



] 



+ Ar [u {n - mx + Ni i t— -t- 6r^ + (n - 1 )N } + 
L 0 l r G r <Sr r n 0 0 rn J 



+ v ( N | | 1 + — - my 1 + N w + I p^ + 

CT |v|r 2 G ; wr 0 yz*0 



9N 

+ q (N - I ) + 2N | | I r I + — 

0 qr zx r I r M 0| 3r 



] 



+ A$ [ (x g WT0T - x q B) cos 0 q cos <}> 0 ] + 

+ A0 [y WTOT - y B) cos 0 - (x W - x B) sin 0 sin A ] + 

u G B 0GB 0 0 J 

+ A6r t (N 6r + N 6rn (T1 0 “ 1 > )u 0 + N |r|6 r U o| r o| ^ + 



3 N 



CRFW 
3 v 



+ An Tn u^r ^ + N u + N | | v / + w rt + N. u 6r 

1 m o o vn o o v v n o o o 6rn o o J 



» , .2 w .... ("0<x> * 2 " 2 > x » 

2 J P( x ) d» > 1 1 * H(X> =======- d * 

/Vq(x) + wjx) 



3N 



CRFW 
3 w 



■ |(-‘ 2 J I 5Khb?) 



v Q (x) w Q (x) 

j 2 2 

/v Q (x) + w Q (x) 



dx 
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3N 



CRFW 

3q 



p ( „ 2 ww ^ f 2 ^ 0 

= 2 ^ J D( ' x) ~ dx ~) U X HU) 7= 



V^(x) w Q (x) 



dx 



/v^(x) + Wg(x) 



3N 



CRFW 
3 r 



- K- t2 ;n*7h*) 



(w 2 (x) + 2 v 2 (x)) 



dx 



/v 2 (x) + w q(x) 
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APPENDIX B 



Table B-l. Nominal design point. 



(FT/SEC) 


v Q (FT/SEC) 


w 0 ( F T / S E C ) 


47 . 330 


2.270 


0.5313 


P 0 (RAD/SEC) 


q Q (RAD /SEC) 


r 0 (RAD/SEC 


-.00003 


0 . 00305 


-.01829 


(DEG) 


<)b 0 (DEG) 


5 r q (DEG) 


- . 7350 


0 . 0000 


2.0000 


F HI 0 (DEG) 


THETA- Cl (DEG) 




-9.4520 


-0.03156 





RPS 

2.369 
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30736E-01 O.V*545E-02 



O.OOOO'JE+OO 0 . OOOOOE *00 O.OOOOOE+OO 0.00000E*00 -.16922E*01 0.9S692E*01 0.72573E-05 0.26399E-03 0.00000E*00 O.OOOOOE^OO O.OOOOOE+OO O.OOOOOE+OO 
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APPENDIX C 
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C-l. 



Nonlinear Simulation of ip - z Compensator for 



Various Commanded Inputs 



(The simulation in this appendix is based on the following table 
of command inputs for heading rate, in units of "tens" of degrees, and 
depth rate, in units of 0.01 * feet/s). 



T i me 


PSIDOT 


Z DOT 


0 . 


- 10.000 


0 . 0 


b 0 * 


- 1 0 .000 


0 . 0 


65 . 


- 1 5 . 000 


0 . 0 


120 . 


- 15. 000 


0 . 0 


1 25 . 


- 15. 000 


0 . 0 5 


1 80 . 


- 15.000 


0 .05 


1 85 . 


0.000 


0 . 0 


240 . 


0 .000 


0 . 0 


245 . 


- 5.000 


0 . 0 


290 . 


- 5 . 000 


0 . 0 


295 . 


- 5 . 000 


- . 1 


350 . 


- 5 . 000 


- . 1 


355 . 


- 10. 000 


0 . 0 


440 . 


- 1 0 . 000 


0 . 0 


445 . 


- 15.000 


- . 1 


6 00 . 


- 1 5 .000 


- . 1 
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Appendix C (Cont . 



CONTROL INPUTS 
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